M Aydinlilar

C Zanni

Fast ray tracing of scale-invariant integral surfaces

: Scale-invariant integral surfaces provide a way to define smooth surfaces from skeletons with prescribed radii defined at their vertices (with linearly interpolated radii along the skeleton edges). The generated surface can be seen as a smoothed version of an infinite union of spheres centered on the skeleton edges. We propose a new rendering pipeline allowing to visualize such surfaces in real-time. We provide comparison to revisited state of the art techniques on a large range of skeleton types for a variety of GPUs. Our method provides improvements for various resolutions on all combination of tested models and GPU hardware (see right graph. More comparisons in the results section).

One of the main difficulties when working with implicit surfaces is to provide efficient visualization. During modeling and final rendering (i.e., processing), visualization may have different requirements in terms of time and approximation. For interactive modeling it is best to have a method requiring no pre-computations and no transformations of the native representation. Provided that short enough computational times can be reached to achieve an acceptable frame rate, ray tracing provides several advantages over meshing. In addition to being trivially parallelizable, ray tracing computations are done in the image space therefore limiting the computations to the parts of the object that are actually visible. This also opens the opportunity to develop output sensitive methods while offering high quality rendering.

When ray tracing implicit surfaces the first challenge is to obtain a robust computation of the closest ray/isosurface intersection along the ray. If we target integral surfaces, such as SCALIS, two additional problems have to be tackled. First, in order to perform field evaluation, one needs to determine the primitives that influence a given location in space. For large skeletons with many primitives, searching through the primitives to find the ones that influence a given location is the main bottleneck. Second, since the SCALIS field evaluations for individual primitives are relatively expensive to compute, a brute force approach requiring a large number of field evaluations is not an option for any interactive application.

In order to overcome these challenges and provide a fast and precise ray-tracing-based integral surface visualization, our strategy combines four main components; 1) We introduce, during rendering, a segmentation of rays in intervals that captures main field variations while generating only a small number of sub-intervals to be processed. 2) We transform field values to a space that allows good quadratic polynomial approximation of the SCALIS field on each of the sub-intervals defined in the previous step along the ray.

3) Fast quadratic roots computation allows to iteratively build polynomial interpolations that rapidly converge toward the real field function. And 4) in order to limit computational overhead, selecting primitives whose supports overlap a given interval along the ray is performed using a dynamic view-dependent acceleration structure.

The last component is built on the fly relying on fast GPU construction of A-buffers which have been previously employed for metaballs rendering [START_REF] Bruckner | Dynamic Visibility-Driven Molecular Surfaces[END_REF] and transparency [START_REF] Thibieroz | Order-Independent Transparency using Per-Pixel Linked Lists[END_REF][START_REF] Maule | A survey of raster-based transparency techniques[END_REF][START_REF] Maule | Memory-Efficient Order-Independent Transparency with Dynamic Fragment Buffer[END_REF]. This results in a method that goes beyond Sherstyuk's fast raytracing algorithm for convolution surfaces [START_REF] Sherstyuk | Fast Ray Tracing of Implicit Surfaces[END_REF], notably allowing a better support of large radius variation while being less sensitive to skeleton tessellation.

Background: Implicit Modeling

In this section, we present required background on skeleton-based implicit modeling and put forth the notations that will be used throughout the paper (for a more detailed introduction to implicit modeling, see [Blo97]).

An implicit surface is defined by a scalar field f and an isovalue c as the set of points p in space satisfying the equation f (p) = c.

Figure 2: Left: a skeleton with prescribed radius on its vertices, Middle: infinite union of spheres defined by linear interpolation of prescribed radius along skeleton edges (this is also a union of sphere-cone), Right: resulting scale-invariant integral surface which can be seen as a smoothing of the sphere-cone union.

Different surfaces can be blended by combining their fields' con-54 tributions. Such combinations are often performed hierarchically 55 in a Blobtree data structure [START_REF] Wyvill | Extending the CSG Tree. Warping, Blending and Boolean Operations in an Implicit Surface Modeling System[END_REF], where each node of the tree 56 represents a composition operator and each leaf is an implicit prim-57 itive.

58

Among implicit surface definitions there exist several conventions on the meaning of the field values. On one hand we have the functional representation F-Rep [START_REF] Pasko | Function representation in geometric modeling: concepts, implementation and applications[END_REF] usually close to a signed-distance field where a volume is defined by f (p) ≤ 0 and where the field f tends toward infinity when moving away from the surface. On the other hand we have density field such as metaballlike surfaces [Bli82; WMW86; NHK*85] defining a volume by f (p) ≥ c (usually with c = 1/2 or c = 1) and where field value tend toward 0 away from the object. Given a decreasing function k called a kernel, the simplest density fields can be defined by applying k to the distance to a point primitive. Such fields can then be combined with composition operators. The most simple smooth blending operator is the sum, used to blend an arbitrary number of input primitives:

f (p) = ∑ i f i (p) (1)
where f i is the field function defined by the i-th primitive. Among kernel families used in practice, we can note the Compact Polynomial kernels that are defined as:

k i,σ (d) =      1 -d σ 2 i 2 i f d < σ , 0 otherwise.
(2)

and the Cauchy kernels which can be defined as:

k i,σ (d) = 1 1 + d σ 2 i 2
(3) They form families of kernels parametrized by a degree parameter i 59 (defining the smoothness of the surface for the compact polynomial segments, triangles or Bézier curves. Such primitives can by orga-1 nized in a skeleton structure in order to facilitate their edition. In 2 this work, we focus on segment primitives.

3

We define a skeleton as a set of line segment primitives S i with 4 prescribed radius information at the endpoints (or vertices -see Fig- 5 ure 2, left). For each point q on the line segment, a radius τ Si (q) can 6 be defined as the linear interpolation of the radii of the two end-7 points. Hence, for each such line segment primitive S i , a sphere-

8
cone can be defined as the infinite union of balls defined by the 9 skeleton points q belonging to the segment with associated radius 10 τ Si (q) (see Figure 2, middle).

11

We present here all the implicit primitives that will be discussed in the remainder of the paper; all these primitives are defined in terms of distance between a point p in space to a point q on the line segment:

d(p, q) = p -q (4)
which, divided with the radius τ Si (q) at the skeleton point q, defines the homothetic distance to a single skeleton point:

h(p, q) = d(p, q) τ Si (q) (5)
or to line segment primitives :

h Si (p) = min q∈Si h(p, q) (6)
Given a kernel k, this allows to define a density function for the line segment S i :

f i (p) = max q∈Si k • h Si (p, q) (7)
with •, the composition operator.

When using such segment primitives with summation blend, bulging appears at skeleton junctions. This problem can be resolved by using integral surfaces, such as a convolution surface [START_REF] Bloomenthal | Convolution surfaces[END_REF] or SCALIS primitive [START_REF] Zanni | SCALe-invariant Integral Surfaces[END_REF]. Both formulations provide an independence from skeleton subdivision when using blending by summation thanks to the additivity property of the integral. Note that both representations are equivalent for constant unit radius. The SCALIS field, main focus of this paper, is defined as:

f i (p) = 1 N k,c Si k • h(p, q) τ Si (q) dq (8)
where the normalization factor N k,c , which depends on the chosen 13 isovalue c and kernel k, is used in order to achieve the prescribed ra- In addition to a direct radius control, this formulation simplifies 18 modeling by providing scale-invariance properties: scaling both the skeleton geometry and the associated radii with a factor s results in 20 a scaling of the isosurface of interest by the same factor -e.g. blend-

21
ing behavior is independent from the scale. This scale property also 22 reduces blurring of details when blended in larger shapes.

23

The method developed in this paper holds for both Equation (7)

24 and (8) with a summation blend. (with linearly interpolated prescribed thickness τ). In this case, the support of an individual primitive (e.g. the volume for which the field value is not zero) is also a sphere-cone that is defined by scaling the two primitive endpoints' radius by a factor σ (see Figure 3 which also describes the convention for the display of the fields).

The Inverse kernel can be seen as a special case of the Cauchy Kernel when the scale parameter σ tends toward 0. In the remainder of the paper, we will discuss the implication of those end-point correctors or their absence. In some cases, such as a truss structure without end-points, all the assumptions related to the presence of end-point correctors can be used.

End-point corrector

Previous work: Raytracing implicit surfaces

Rendering implicit surfaces is a subtle trade-off between efficiency and accuracy and has triggered a lot of research. In order to visualize an implicit surface with the ray tracing approach, one has to find the first parameter along a given ray:

r(x) = o + x d (9)
submitted to COMPUTER GRAPHICS Forum (2/2021).

for which the field value evaluates to the isovalue of interest c, e.g. finding x such that :

f • r(x) = c
or equivalently :

f • r(x) -c = 0 (10)
Henceforth, we only consider Equation (10) and discuss the detection of zero crossings along the ray.

A first solution to do so it to apply the ray marching approach to isolate the root by advancing with a constant step size until a sign change is detected in f • rc. Then, applying another method, such as a constrained Newton method, to localize the root more precisely. The behavior of such method is highly correlated to the choosen stepsize: for small steps the method becomes computationally expensive, for larger step size, some isosurface crossing can be missed.

Available techniques tend to differ depending on the properties of the implicit surface to be ray-traced, for instance density field versus distance field.

Polynomial definitions/approximations Several ray tracing algorithms rely on polynomial root finding algorithms to locate the isosurface along the ray, either because the field itself is defined by polynomials or because a polynomial approximation of the field is built along the ray.

For point primitives defined with the compact polynomial kernel (Equation (2)), the field function along the ray is a piecewise polynomial. For small degree (up to degree 4), this property can be used Choosing an adequate level of subdivision is problematic: on one hand, small number of intervals can result in poor approximation depending on both the viewpoint and the skeleton tessellation even more so in the context of SCALIS as it becomes more difficult to capture the maximal field contribution of a given primitive. On the other hand higher number of intervals results in higher computational time due to additional field evaluation. Furthermore, when performing computation on GPU, the incoherence of the generated subintervals require additional data management to avoid the computation of additional field values.

An alternative to interval subdivision is to rely on a higher degree approximation defined from additional sampling positions [She99;

in the gradient and the sign of the field function. Such approach The approximations are then used to compute an approximation of the isosurface position. To increase the approximation precision, supports are uniformly subdivided before generating the final segmentation which can result in intervals of highly incoherent size. Middle: our ray segmentation only relies on estimated maximal influence of individual primitives to define an initial segmentation with a small number of subintervals. We then rely on dynamic subdivision in order to increase precision of approximation wherever required. Right: Hermite data defined at interval's end-points and associated polynomial interpolation. visualization algorithm [START_REF] Sigg | GPU-based Ray-casting of Quadratic Surfaces[END_REF]. We extend this approach for better management of segment primitive with varying radius (primitives whose support are also sphere cones).

Our approach

Our main shape to be ray-traced is the SCALIS field (see Equation (8)) defined by a skeleton consisting of line segment primitives. Since SCALIS field evaluations are expensive to compute, brute force ray-marching approaches are not feasible. We therefore propose to rely on local field interpolations in order to reduce the number of field evaluations required for calculating the ray-surface intersections. To achieve high fidelity interpolations with lower degree polynomials, we propose to subdivide the ray into intervals that capture the main field variations along the ray. A different interpolation is used and refined in each interval, locally producing an accurate approximation. The overall idea and the processing of an example ray can be seen on Figure 5.

Our approach relies on a bijective mapping of the SCALIS field to an approximate smooth homothetic distance field (squared) which has the exact same isosurface of interest as the SCALIS field.

All processing is done on this normalized field. This new field exhibits similar variations as the homothetic distance fields (squared)

to individual primitives h Si 2 (see Equation ( 6)).

We use the correspondence between the normalized field and the h Si 2 fields in order to define our ray-subdivision strategy: the intervals are defined by cutting the ray at the minima of h Si 2 for each individual primitive S i (purple curves in Figure 5 (a)). We expect such cuts to limit the number of oscillations of the normalized field in a given subinterval, a property required to perform field analysis through low degree polynomial interpolation. As can be observed in Figure 5 (a), the normalized field has this desirable property on each sub-intervals. By cutting the ray at the minima of h Si 2 , the likelihood of finding a point within the implicit surface is increased, such points allows to isolate a root by detecting a sign change. While the local minima of the homothetic distance field along the ray provides only an approximation to the local minima of the SCALIS field, this approach behaves well in practice (see section 8).

Our ray-tracing algorithm generates and processes intervals on the fly in depth-order until an isovalue crossing is detected. This corresponds to the main loop (A) in the diagram in Figure 5. We process a subinterval by iteratively refining a local polynomial interpolation. The interpolation is first initialized using the Hermite data sampled at interval end-points. Then, it is refined by reducing the interval extent, cutting it at the root of the polynomial which interpolates the field on the reduced interval (inner loop (B) in the diagram and graph of Figure 5). Provided that the polynomial interpolation has a root within the sub-interval, the root of the successive interpolations will converge toward the real isosurface as the interval gets iteratively smaller. We use smooth piecewise quadratic polynomials to obtain fast and stable root computations.

We first discuss the mapping from density field to homothetic squared distance field in section 4. Then, in section 5 we describe our ray subdivision strategy and in section 6 we present the approach used to perform the processing of a ray's subinterval. Fi- nally, we describe in section 7 our GPU implementation relying on a dynamic data structure.

Approximated squared homothetic distance

Scale-invariant integral surfaces are defined from the homothetic distance to skeleton points as described in section 1.1. For most of the kernels used in practice, it is actually defined from a squared homothetic distance. We use this property to introduce a mapping that allows to compute more precise quadratic interpolations. We first present the specific configurations where a quadratic polynomial can be fitted exactly to the mapped values, then discuss the general case.

Remapping SCALIS field values

In previous works on the integral surfaces, the study of infinite primitives has proven its usefulness for thickness control around skeletons [START_REF] Zanni | SCALe-invariant Integral Surfaces[END_REF] and blending control [START_REF] Zanni | Nary implicit blends with topology control[END_REF]. We use it to analyze the homothetic distance to the isosurface and provide better polynomial interpolation. We first introduce a field mapping for infinite line primitives, then generalize it for arbitrary skeletons.

Infinite line primitive Let us consider an infinite line primitive with a constant prescribed radius τ 0 . For the main kernel families used to define skeleton-based implicit surfaces, the SCALIS field for such primitive in isolation can be defined as a function of distance d to the line by using a kernel k of the same family with a different degree. Hence, the field can be defined as :

f line,τ0 (d) = f line,1 d τ 0 = λ k d τ 0 ( 11
)
where λ = c/ k(1). For the Compact Polynomial kernel of order i,

we have k = k i+1,σ and for the Cauchy kernel, we have k = k i-1,σ .

20

We generalize the approach of [START_REF] Bruckner | Dynamic Visibility-Driven Molecular Surfaces[END_REF] from metaballs to full skeletons with varying radius. The homothetic squared distance (instead of a Euclidian distance) to the line primitive (instead of point primitive) can be computed at any given location in space from the field value by applying the inverse of the function defined in Equation (11) followed by squaring the result :

g( f ) = f line,1 -1 ( f ) 2 = k-1 f λ 2 = d τ 0 2 (12)
For the Compact Polynomial kernel, this develops into:

g( f ) = d τ 0 2 = σ 2 ((1 -(1 -1 σ 2 )( f c ) 2 i+1 )) if f > 0 , 1 otherwise.
(13) while for the Cauchy kernel, this develops into:

g( f ) = (1 + σ 2 ) c f 2 i-1 -σ 2 (14)
with the scale parameter σ and the order i of the kernel as defined in Equation (2,3). Note that in the case of radius τ 0 = 1, this is also the isosurface of interest can be computed:

g( f ) -1 (15)
In the remainder of the text, we call the field resulting from the application of Equation ( 15) to the SCALIS field the normalized field.

For a line with constant radius in isolation (or a long enough segment for the compact polynomial kernel), it is important to note that for a given ray, the squared distance from the ray to the line is a degree two polynomial. Hence, for such a configuration, the squared homothetic signed distance can be exactly interpolated by a quadratic polynomial defined by Hermite data sampled at an arbitrary position along the ray (i.e. by evaluating g • f • r -1 and its derivative at the ray parameter t).

Other special configurations For a few other configurations, the Generalization Interestingly, Equation ( 13) can be used directly on the SCALIS field generated from a more general skeleton in order to compute a smooth approximation of the homothetic squared distance to the surface. Indeed, since g is a strictly decreasing function on R + , hence injective, the normalized field leaves all the level set unchanged (in terms of geometry), including the one corresponding to the isosurface of interest (which now corresponds to an isovalue of 0 instead of c). And, outside of the blending area (e.g. branching) and the neighborhood of skeleton end-points, the field behavior will tend toward one of an infinite line before the application of g, hence toward a squared homothetic distance after application of g.

function g • f • r -1 is
We designed the normalized field to have good properties for our algorithm: it does not create new oscillations in the field and is less dependent on the kernel used. The mapping has the beneficial property of removing the inflexion points of f • rc that are due to the Cauchy and Compact polynomial kernel shape (see Figure 6).

Ray subdivision

In order to simplify the processing of a given ray, we want to define intervals with limited number of oscillations (typically either one or two local minima of g • f • r -1 exists in the defined intervals).

This provides better interpolations in configurations where the normalized field along the ray is not exactly a quadratic function.

We further use the fact that the normalized field can be seen as a smooth approximation of the homothetic distance (squared) to the surface min i (h 2 Si ) -1 in order to simplify the processing of a given ray. A high correlation between these values can also be observed, and is even more noticeable after mapping f • r to approximate homothetic squared distances (blue curve).

Correlation with h 2

Si -1 50

In order to ease the analysis of the field, we take inspiration from 65

The aforementioned evidence hints that a good sampling strategy 66 is achievable by computing the arguments of the minima of the 67 homothetic distance to the individual primitives h Si • r.

Argminimum of homothetic distance along a ray

1

In order to generate our ray subdivision on the fly, we need to find efficiently the argument of the minimum of the homothetic distance from a segment primitive S i with linearly varying radius (Equation (6)) to the ray r, i.e. we are trying to solve:

argmin t h Si • r(t) (16)
Let us define the segment primitive S i , parametrized by:

p(s) = p 0 + s u and τ(s) = τ 0 + s ∆τ (17)
with u a unit vector and the parameter s in the range [0, L], L being the length of the segment. For a skeleton point (p(s), τ(s)) in isolation, the Euclidean distance and the homothetic distance have the same argument of the minimum along the ray, the minima only differing by a factor 1/τ(s). For the Euclidian distance, the minima (squared) have a closed form expression :

dr 2 (p(s)) = (m + su) -(m + su) T d d 2 (18)
with m = p 0 -o, where o is the ray origin and d is the ray direction.

2 Therefore, in order to solve Equation ( 16), we can invert the order in which minima are computed and we can instead study:

argmin s∈[0,L] d 2 r (p(s)) τ(s) 2 (19)
This gives the skeleton parameter s min for which the minimum end-points. This problem is related to the problem of radius shrinkage at skeleton end-point observed in [START_REF] Zanni | SCALe-invariant Integral Surfaces[END_REF]. When applying the end-point correctors, we force an additional sampling position near skeleton end-points, increasing the correlation between the argument of the minima of h Si • r for individual primitives and the ones of g • f • r -1. Note that in the absence of correctors, a heuristic presented in section 6.2 can compensate this shortcoming.

Ray processing

The ray subdivision strategy presented in the previous section not only provides intervals with limited number of oscillations, it also limits the initial number of intervals along a given ray. Smaller number of intervals result in smaller number of compulsory field evaluations while processing a given ray.

We also subdivide the ray when switching from an area not overlapped by any primitive support to an area overlapped by at least one primitive support (and vice-versa). This way, we can avoid processing large empty areas. This is easily achieved with the data structure presented in section 7.

We first present the quadratic polynomial interpolation used to define local field interpolation within a given subinterval, then we present the processing of a given subinterval.

Quadratic polynomial interpolation

Working with quadratic polynomial interpolation has the advantage of much more efficient, simpler and numerically stable root computation. As discussed in section 4.1, thanks to the field normalization, it also provides exact interpolation, in specific cases, independently of the kernel used. However, some configurations of Hermite data at interval's end-points cannot be realized by a quadratic interpolation (four constraints for only three degrees of freedom, see Figure 7). In such incompatible configuration, a possible approach would be to rely on Hermite-Birkhoff data, ignoring one of the two derivatives at interval end-points (see green curve in Figure 7). Instead, we propose an alternative solution that allows to respect the Hermite data at both end-points by relying on two piecewise quadratic polynomials, and without computing any additional field values. For simplicity we cut the initial interval in two sub-parts of equal size. Our interpolant is defined as a one dimensional Bézier curve on each part. For such curves, the control polygons are constrained such that the control points are equally spaced in abscissa. The ordinates of control points are chosen to verify Hermite data and interpolation smoothness. From the six control points Q i=0..5 (three for each of the Bezier curves), the two first control points Q 0,1 of the first interval (respectively, last points Q 4,5 for the second intervals) are constrained by Hermite data. Remaining control points on both sub-intervals should have the same value Q 2 = Q 3 to ensure continuity and they should lie on the line joining the two adjacent control points Q 1 and Q 4 to ensure smoothness, which leads to the following ordinates:

1 2 (h 1 + h 2 ) + L 8 (h 1 -h 2 ), (20)
where (h i , h i ) is the Hermite data at each end of the interval. Note 1 that this value is also equal to the value returned by the cubic poly-2 nomial interpolation at the middle of the range (see Figure 7) as For the Cauchy kernel, we obtain much better interpolations with the normalized field for small numbers of interval's subdivision, as such kernel is badly suited for low degree polynomial interpolation on larger intervals.

Processing of an interval

We generate intervals iteratively and process them on the fly. For each interval we apply a routine in order to check if the isosurface is crossed. If so, we return the position of the isosurface along the ray, else we move to the next interval (see section 7).

Let us define [t begin ,t end ] as a sub-interval to process. We use the quadratic interpolation of Hermite data presented in Section 6.1 in order to estimate the isosurface position if it exists. The estimated position troot is defined as the first root of the interpolation defined from Hermite data evaluated at t begin and t end .

The routine iteratively reduces the size of the interval in order to define Hermite interpolations of increasing precision (see Fig- We rely on rational Bézier curves until a root is guaranteed. Such 5 curves allow to have an interpolation that is arbitrarily closer to 6 the control polygon hence increasing the chance to detect an iso-7 value crossing (see Figure 9). Note that with rational functions, primitives supports are stored per ray, then processing the ray. We 25 rely on the stored data to serve as an acceleration data structure to efficiently retrieve the primitives whose supports overlap the given interval of the ray. When using the Cauchy kernel, we clip the support of a given primitive by a sphere-cone whose scaling σ clipping is computed to guarantee a maximal error ε clipping . This allows us to use the same approach for both kernel families.

The first step relies on a fast GPU construction of Abuffers [Thi11; MCTB11; MCTB12]. A geometry shader is used to generate a quad per skeleton segment such that the quad covers the projected segment primitive support. As the support of a given primitive is defined by intersecting some quadrics (a sphere-cone), we compute the quad using [START_REF] Sigg | GPU-based Ray-casting of Quadratic Surfaces[END_REF]. We adapt the proposed approach to compute segment-aligned quads to limit the number of generated fragments for which ray/sphere-cone intersections are computed (in particular for primitives whose maximal radii are comparatively small to the segment length). The ray-supports intersections are computed in a fragment shader. Entry and exit points in the sphere-cones, as well as the associated primitive id, are inserted in the A-buffer linked-list (see Figure 10) and sorted on the fly according to the depth of the entry [START_REF] Lefebvre | Per-Pixel Lists for Single Pass A-Buffer[END_REF]. Note that contrary to [START_REF] Bruckner | Dynamic Visibility-Driven Molecular Surfaces[END_REF], the use of a compact support kernel avoids the need to cut-off primitives influences in this first step.

Of course, any variant of A-buffer creation could be used to accelerate this first step, for instance postponing the sorting of fragments until the ray processing step [START_REF] Bruckner | Dynamic Visibility-Driven Molecular Surfaces[END_REF] or using advanced memory management [START_REF] Schollmeyer | Direct Isosurface Ray Casting of NURBS-Based Isogeometric Analysis[END_REF].

Once the first step of the processing is done, the linked-lists of the A-Buffer define intervals with a fixed set of segment primitives. We now face a choice. Subdividing the ray based on those intervals would minimize the maximum number of segment primitive to be stored per fragment shader call. However, it would also increase the number of field computations required. Increasing the interval size too much might require more memory to store all segment primitives whose supports overlap each intervals. On the opposite, subdividing the ray only when switching from an area not Figure 10: A-buffer linked-list associated to a given pixel after rasterization of all skeleton's primitives. Note that entry depth z in and exit depth zout in the support of a given primitive are stored in the same linked-list node. This linked-list is the main entry to the ray processing algorithm and allows to retrieve efficiently the primitives whose supports overlap a given ray's subinterval. overlapped by any primitive support to an area overlapped by at least one primitive support (and vice-versa) would require storage of a large number of segment primitives at once depending on the ray/skeleton configuration, hence increasing memory usage of each individual thread.

We choose to combine this second subdivision strategy to our contribution presented in section 3. This presents a nice trade-off: it limits the number of segment primitives to be stored at once while generating long enough intervals, hence limiting the number of compulsory field evaluation. The fragment shader that renders the view generates and processes the intervals on-the-fly using algorithm 1 until an intersection with the isosurface is found. The resulting linked-list processing loop is given in Algorithm 2.

Due to possible numerical instability in the normalized field derivative computation (the latter being non continuous when the field switch from null to non null), we apply one step of sphere tracing to avoid those areas, note that this require no field evaluation as the field is null at the support limit (see comment in Algorithm 2).

Toward output sensitivity Similarly to Bruckner's approach to achieve output sensitivity [START_REF] Bruckner | Dynamic Visibility-Driven Molecular Surfaces[END_REF], it is possible to render an occluding depth buffer in order to limit the number of primitives that will be registered in the linked-lists. The depth buffer is initialized with depth values corresponding to positions that are guaranteed to be in the volume. Then all primitives whose support entries are beyond the occluding depth can be safely ignored. For efficiency those depth values need to be computed in a first render path on a per primitive basis. Due to the integral formulation of the field, the minimal radius around primitives does not necessarily correspond to the prescribed radius (e.g. typically for bent skeletons). In order to overcome this difficulty, we rely on two assumptions: usage of the skeleton end-point correctors and bounded radius variations: none of the two spheres at a primitive end-point completely encompass the other -e.g. each vertex should have an influence on the geometry of the union of sphere. Note that we do not render occluders for primitives corresponding to end-point correctors. In this context, thanks to the scale invariant property of the SCALIS representation, it is possible to analyze the worst skeleton configuration to derive a scaling factor for the sphere cones' radii such that the implicit surface is guaranteed to enclose the modified primitives. the blue curve corresponds to constant radius). In section 8, we discuss rendering times with and without this optimization. In absence of end-point correctors, an alternative strategy would be to rely on the field behavior of primitives in isolation. Indeed, for long enough primitives with constant radius, there is a portion of the segment primitives for which the prescribed radius is guaranteed. While we 57 do not use this property in our implementation, it could be used to 58 improve the occluders for some specific skeleton configurations. Note that in practice, it is is also possible to blend different 16 skeletons with a blending sharper than the summation blending, 17 e.g. by relying on Ricci's blending [START_REF] Ricci | Constructive geometry for computer graphics[END_REF] (see Figure 14). Indeed, ing to obtain visually acceptable results on primitives with constant 29 radius. Our implementation of the sphere-tracing routine also uses a ray subdivision strategy in order to discard space outside of any primitive support as well as to limit the maximal number of segments to be stored and processed at the same time during the ray processing loop. On each subinterval the number of sphere-tracing steps is bounded (typically 256 in our implementation).

For the compact polynomial kernel, all renderings are done at a resolution of 2048 × 2048 and rendering times are provided in Table 1. On all tested examples, we observe improvements of the global rendering time with our method in comparison to the most competitive one (which is dependent on the model). Rendering time reduction ranges from -34.2% to -48.7% on the Intel card, from -51.5% to -78.0% on the nVidia Quadro card and from -35.3% to -72.1% on the nVidia RTX card. In addition, both sphere-tracing and Sherstyuk's fast ray tracing can present visual artifacts (see Figure 17, note that those artifacts can be reduced at the expense of additional computation time). It is also important to note that in presence of a large amount of blending -for instance if a large number of segments connect in one vertex (such as in Figure 12) -the sphere-tracing algorithms would require usage of under-relaxation to properly find the surface.

We also have tested our method at higher resolutions where we observe similar improvements (see Figure 20). In order to perform a stress case on the nVidia RTX, we rendered a foam structure consisting of 526k segments using the compact polynomial kernel, running at 15 frames per second (see Figure 12 In each direction, we sample the ray origins uniformly in the asso-

24

We observe a reduction for all statistics when using our method (see Table 2), including when primitives present small radius variations (e.g. the radius is either approximately constant over the full skeleton or is approximately constant by part on the skeleton). Larger reductions are observed in presence of varying radii. However, the main benefit is achieved along grazing rays that present the maximal number of steps (a well-known problem of the sphere-tracing algorithm -see Figures 21 and18).

When using distance primitives instead of SCALIS primitives, this limitation of the sphere-tracing algorithm can be alleviated by the use of directional Lipschitz bounds (e.g. the segment-tracing algorithms). In this context, our method still provides reduction of the maximal number of steps required along grazing rays (see Table 3).

In order to validate the robustness of our method, we have also measured errors along rays. In practice, less than 10 -3 % of the rays are missing an isosurface crossing, all of them are grazing rays (see Figure 18 and 21). Those few errors are due to Hermite data configurations for which control polygons of the interpolant do not cross the isovalue 0. This limits the efficiency of our rational interpolation heuristic. In order to mitigate this shortcoming, it could be interesting to investigate 2D quadric Bezier curves which would allow to increase the amplitude of the control polygon by moving the abscissa of the control points.

Limitation Our algorithm and its implementation present a few limitations. First, the heap used to store segment primitives during the ray processing loop should be large enough to accommodate all primitives whose supports overlap the initial intervals defined by the argument of the minimum of homothetic distance. Observe that this limitation could be mitigated by adding an additional break condition in Algorithm 2 based on the current heap size. Secondly, output invariance could be improved by combining our current strategy with occluders generated from primitives considered in isolation. Finally, the A-buffer can become expensive to build for finely tessellated curves (due to overlaps between primtive supports). Two main directions could be investigated : curve primitives [START_REF] Fuentes Suárez | Anisotropic convolution surfaces[END_REF] and curved simplification based on local radius and kernel scale. Similarly, finely tessellated skeleton also increase the number of argument of the minimum of homothetic distance along the ray.

Conclusion and Future work

We have introduced a new rendering algorithm for scale-invariant integral surfaces. Our approach presents no visible artifacts while decreasing computation time by up to 70% in comparison to revisited previous work. Our key contributions are 1) the usage of a new initial sampling strategy of the rays based on homothetic distance to segment primitives, 2) the introduction of a new field mapping combined to quadratic polynomial interpolation, and 3) an efficient GPU implementation relying on a A-buffer data structure. We believe this technique could help in further study of a larger range of integral surfaces.

Among future directions of research, we can mention management of triangle primitives (these share the main properties used in our method) and transparent rendering (relying on depth slabing to submitted to COMPUTER GRAPHICS Forum (2/2021).

Object

Intel [START_REF] Peytavie | Arches: a Framework for Modeling Complex Terrains[END_REF] and our method on blob primitives.

limit A-buffer creation cost and memory usage). Finally, we believe 1 that a combination of an A-buffer acceleration structure and/or a

14 dius

 14 around the skeleton. Visually, the resulting isosurface appears 15 as a smoothing of the sphere-cones associated to the prescribed 16 radii (see Figure 2, right).

Figure 3 :

 3 Figure 3: Leftmost: line segment primitives with prescribed radius at vertices. Left: associated sphere-cones. Right: Primitive supports when using a compact kernel. The supports are also spherecones and their radii are proportional to the radii τ Si prescribed on the skeleton. Rightmost: a slice of the SCALIS scalar field, the points that are inside of the volume defined by f (p) > c are depicted in red, the points outside are depicted in blue. The points outside of the supports of all primitives have a null field value and are depicted in white.

  Due to the definition of the field as an integral, see Equation (8), the prescribed radius is not achieved in the tangential direction near skeleton end-points. In [ZBQC13], a simple correction scheme is presented to overcome this drawback. It consists in adding a segment primitive with constant radius, the length of which is proportional to the prescribed radius at the end vertices of the skeleton as well as on a constant depending on the kernel used (see Figure inset).

  to compute a closed-form expression of the root [GPP*10]. For higher degrees, the Bézier clipping algorithm can be used[START_REF] Nishita | A Method for Displaying Metaballs by using Bezier Clipping[END_REF][START_REF] Kanamori | GPU-based Fast Ray Casting for a Large Number of Metaballs[END_REF] to iteratively converge toward the first root of the polynomial (or to rapidly reject root-free intervals).For convolution surfaces, polynomial approximations of the field value can be built along the ray[START_REF] Sherstyuk | Fast Ray Tracing of Implicit Surfaces[END_REF]. During rendering, each primitive can be approximated by one or several polynomials on the interval defined by the intersection between the ray and the primitive support. The interval is uniformly subdivided and polynomial approximations are calculated from Hermite data (field value and derivative) sampled at subinterval boundaries (see Figure4). The roots of resulting cubic polynomial approximations can be computed analytically. When primitives are combined with the summation operator, new polynomials are defined by summation of polynomials on the intersection of primitive subintervals.

71Figure 4 :

 4 Figure 4: Left: In [She99], boundary of primitive supports are used to subdivide the ray in order to compute local polynomial approximations of field values.The approximations are then used to compute an approximation of the isosurface position. To increase the approximation precision, supports are uniformly subdivided before generating the final segmentation which can result in intervals of highly incoherent size. Middle: our ray segmentation only relies on estimated maximal influence of individual primitives to define an initial segmentation with a small number of subintervals. We then rely on dynamic subdivision in order to increase precision of approximation wherever required. Right: Hermite data defined at interval's end-points and associated polynomial interpolation.

28

  Interval/Affine arithmetics For arbitrary function-based field 29 definition, a robust approach to ray tracing is to rely on iterative 30 interval refinement and interval arithmetics [Mit90; CHMS00] to 31 compute a bound on field variation on a given interval. Each time 32 the estimated bound contains the isosurface c the interval is subdi-33 vided into two. Affine arithmetics [FPC10] can be used to obtain 34 tighter bounds. 35 In [Kee20], one of the main drawback of this category of meth-36 ods is tackled (albeit on an alternative direct rendering approach): 37 efficient computation on large expressions. It relies on a represen-38 tation of the expression well adapted to GPU evaluation as well 39 as on-the-fly simplification of the expression. This method targets 40 general implicit surfaces, the author hints at reduced efficiency of 41 the interval evaluations in presence of shapes with many smooth 42 blends. Our focus, skeletons with large number of primitives, falls 43 in this category. 44 Acceleration data structures For large skeletons, computational times are also highly correlated to the time required to iterate over all the primitives that influence a given point along the ray (e.g. 47 during field evaluation). 48 In [GPP*10] a bounding volume hierarchy (adapted for the man-49 agement of the blending operation) was used for efficient evalua-50 tion. More recently, [Bru19] has relied on a dynamic data structure 51 that is efficiently rebuilt at each frame by using the GPU rendering 52 pipeline: for each ray, a linked list of intersection points with sup-53 ports of each primitive is computed relying on a GPU based quadric 54

Figure 5 :

 5 Figure 5: General pipeline of our method. Top Left: A slice of the SCALIS field consisting of three primitives and a ray (green) to be processed. Below it, the SCALIS field variation along the ray f • r (orange) and the normalized field values (in blue) are given. We are interested in locating the zero-crossing (root of the Equation (10)) marked with the red disk. In the middle the two main loops of the algorithm are denoted as (A) and (B) and on the right these are shown respectively for the given ray: dividing the ray into intervals (a) and iterative root refinement by polynomial interpolation (b). The normalized field exhibits limited variations on sub-intervals defined by cutting the ray at local minima of the homothetic (squared) distance to individual primitives (purple curves in (a)).

22the

  Euclidean squared distance to the line.Similarly, as g(c) = 1, a signed squared homothetic distance to submitted to COMPUTER GRAPHICS Forum (2/2021).

  also a quadratic polynomial, hence it can be exactly interpolated by a quadratic polynomial defined by any Hermite data sampled along the ray. Such configurations are: any rays for two equal line primitives, all rays belonging to the bisecting plane between two parallel line primitives, the three rays that correspond to the equidistant lines for two crossing line primitives and the two rays corresponding to equidistant lines for two arbitrary line primitives. Indeed in all those cases, the global field is defined by 2 f line,τ0 (d), hence a constant can be factored out from Equation (13).

Figure 6 :

 6 Figure 6: A sample ray on the skeleton defined in Figure 2. Left: Cauchy Kernel. Right: Compact Polynomial Kernel. The field variation along the ray f • r (orange curve) exhibits limited variations on sub-intervals defined by cutting the ray at local minima of homothetic (squared) distance to individual primitives (purple curves) .A high correlation between these values can also be observed, and is even more noticeable after mapping f • r to approximate homothetic squared distances (blue curve).

51the

  Linderberg principle of absence of local extrema creation by 52 convolution with a specific family of kernel [Lin91]. While we are 53 not exactly in this context (Cauchy and Gaussian verify exactly the 54 condition, compact kernel is only close to it, SCALIS is not a con-55 volution when using varying radii), we experimentally observe a 56 similar behavior. Furthermore, the mapping g is decreasing, it can-57 not create new local extrema, which allows to do a similar obser-58 vation on the normalized field (see Figure 5 and 6). Similarly, we 59 expect the local normalized field behavior to directly correlate to 60 homothetic squared distance h 2 Si to individual primitives including 61 in presence of varying radius. Indeed, it is defined by blending (i.e. 62 smoothing) contributions of all skeleton-points based on their prox-63 imity (i.e. with sharper kernels the correlation between g • f • r -1 and h 2 Si • r -1 increases).

36

  is reached. The ray parameter t min can easily be derived from it 4 (e.g. by computing the orthogonal projection of the skeleton point 5 p(s min ) onto the ray). We therefore have to study a rational function whose poles are not in the range of interest (i.e. in the range [0, L]), thus the minimum of the function is either reached at the boundary of the interval [0, L] or where the derivative cancels. By the cancellation of the derivative, we get a linear equation in the form bx + c = 0 with: b = τ 0 L 2 (1 -(u T d) 2 ) -∆τLu T (m -(m T d)d) and c = ∆τ m -(m T d)d 2 + τ 0 Lu T (m -(m T d)d) Note that Equation (16) could be solved directly. Recall that the sphere-8 union associated to a given primitive 9 includes a cone section. Depending on 10 the ray/cone configuration, it amounts to 11 finding the smallest scaling factor to ap-12 ply to the segment's radii such that there 13 exists a single intersection between the ray and the scaled sphere-14 cone, see Figure inset. We find that our approach requires less computations and is easier to compute in a numerically stable way. 16 End-point correctors The minima of h Si • r are highly correlated 17 to the maximal field contributions for a given primitive. However, 18 due to the integral nature of the SCALIS field definition, field vari-19 ations diverge from this behavior near primitive end-points. While 20 this is not problematic in general (due to blending by summation 21 of adjacent primitives), this could become a problem near skeleton

Figure 7 :

 7 Figure 7: Two Hermite data configurations and associated polynomial interpolations: Hermite cubic interpolation in blue, quadratic interpolation by part in red and quadratic Hermite-Birkhoff interpolation in green.

Figure 8 :

 8 Fig.11(3)

36

  figurations, we experimentally observe a similar convergence rate 8

  ure 5). By doing so, the estimate troot converges toward the exact isosurface location.The interval size is reduced by moving one of its endpoints to the estimate troot depending on the Hermite data evaluated in troot .If no clear choice is possible (e.g. no sign change is detected in normalized field value and neither of the sub-intervals [t begin ,tt root ] and [troot ,t end ] has a constant sign for their directional derivative, the first sub-interval is chosen and the second interval is stored to allow one step of backtracking.The iteration continues until the presence of the isosurface can be rejected or the isosurface has been located within a given error threshold. The structure of the interval processing loop is given in Algorithm 1.Heuristic/Oracle On intervals presenting large radius variations (or in absence of end-point correctors) the polynomial interpolation might overestimate normalized field values g• f •r-1, opening the possibility of missing the isosurface. One solution to this problem would be to uniformly subdivide intervals in order to obtain better initial polynomial interpolations. However this comes at a large increase in the number of intervals to be processed. In order to avoid

  this, our algorithm uses a heuristic to limit the need for extra sub-1 division. We modify the oracle estimating the existence of a root 2 (and its position if it exists) until an isosurface crossing has been 3 isolated (i.e. until a sign change exist between interval end-points).

  4

8

  root computation remains the same as long as the pole is not in the 9 range of interest. In our implementation, a weight equal to 3 on the 10 middle control point of the interpolation curve was sufficient on all 11 tested examples. 12 With the modified oracle, even without any interval subdivision, 13 failure cases are rare (see section 8), especially compared with the 14 behavior of the sphere-tracing algorithm that requires arbitrarily 15 large number of steps therefore creating holes in the surface, when 16 working with limited number of steps for rapid rendering (see again 17 section 8).

18 7.

 18 Dynamic data structure : efficient GPU implementation 19 When evaluating the SCALIS field on large skeletons with com-20 pact kernels, the computation time is mostly driven by the selection 21 of primitives whose supports overlap the evaluation point. Simi-22 larly to [Bru19], our GPU implementation works in two steps, first, 23 building a dynamic data structure where entry and exit points in

  24

  submitted to COMPUTER GRAPHICS Forum (2/2021).

  Scaling factor Computation of this scaling factor is discussed in Appendix A. The quality of the minimal volume is related to the kernel scale parameter, the larger the scale the smaller the guaranteed volume is. For instance, for typical Compact polynomial kernel scale σ ranging from 1.5 to 2, the scaling factor range from ≈ 0.62 to ≈ 0.35 (see Figure inset -the purple curve correspond to our bounded radius variation). Similarly, for more constrained radius variations, better bounds can be computed (see Figure inset

  ,

60

  Our method was implemented in C++ using the OpenGL library 61 and shaders were programmed in GLSL. We have tested our im-62 plementation on a large range of objects including artistic shapes 63 (see Figures 1, 2, 13, 15), truss structures (see Figure 16), branch-64 ing structures (see Figure 11), procedural foam structures (see Figand high-end graphic cards (namely an Intel UHD Graphics 630, 67 an nVidia Quadro P1000 and an nVidia GeForce 2080 RTX). Ex-68 amples have been chosen to present a large variety of skeletons 69 in terms of number of segments, local density of primitives and 70 whether or not the radii change. In order to assess the efficiency

Figure 11 :

 11 Figure 11: Left: Procedural trees. Right: Random skeletons with constant and varying radii used to test the resilience of our algorithm.

Figure 12 :

 12 Figure 12: Procedural foams whose skeletons are defined as edges of a voronoi diagram of points.

Figure 13 :

 13 Figure 13: Real-time rendering allows users to efficiently explore kernel scale parameters for a given skeleton.

Figure 14 :

 14 Figure 14: Left: Model generated using a Cauchy kernel and Ricci blendings of subcomponents (main body, eyes, tongue, legs and branch), Middle: Comparison with summation-only blending, Right: Comparison with compact polynomial kernel. .

Figure 15 :

 15 Figure 15: Skeleton-based free form modeling.

Figure 16 :

 16 Figure 16: Truss structure.

  15

18 such

 18 blending keeps the correlation of the field with the homothetic 19 squared distance to individual primitives by defining a blending be-20 havior parametrized between a summation blend and a maximum 21 blend. 22 GPU rendering time We record the average rendering time for a 23 set of viewpoints around our test objects. In order to provide fair 24 comparisons, the GPU implementation of the sphere-tracing and 25 the one of the Shertyuk ray tracing both rely on the A-buffer. In our 26 implementation of Sherstyuk's method, we use four subdivisions 27 per primitive, which correspond to the coarsest subdivision allow-28

  , right). Finally, we provide individual timing for the A-buffer creation step (with and without occluders) and the ray processing steps. As expected, usage of occluders can have a large impact in terms of rendering time (as well as memory consumption) depending on the nature of the skeleton.For the Cauchy kernel, we mainly tested our method on the nVidia RTX card with a lower resolution of 1536 × 1536 in order to account for the higher number of primitive overlaps generated by the larger clipping sphere-cone (e.g. σ clipping = 4 for σ = 0.3 and ε clipping = 0.005), otherwise this would lead to an A-buffer overflow. The use of a kernel with a larger footprint results in an increase of computational time. With our method, timings range from 7 to 54 milliseconds (excluding Fig.12(Middle)). With respect to sphere tracing, this corresponds to a rendering time reduction ranging from -46.2% to -81.3%. We observed similar improvements on the two other GPUs when run at a resolution of 1024 × 1024.

Figure 17 :

 17 Figure17: Left: our method presents minimal errors in isosurface intersections, Middle: when using sphere tracing with a number of steps limited to a few hundreds (in order to limit runtime cost) holes can appear in the surface,Right: Sherstyuk fast ray tracing can produce large deformation of the surface depending on the viewpoint.

Figure 18 :3

 18 Figure 18: Number of field value computations used during the processing of a ray. The red pixels correspond to isosurface crossing missed by our method. Left: Sphere-tracing Right: Our method.

9

  Statistics In order to compute per ray statistics, we have launched 10 sets of rays in six directions (three main axis with both posi-11 tive and negative orientations). For each model, we have com-12 puted the axis-aligned bounding-box from sphere-cone supports.

  13

14

  ciated bounding-box's face such that a total of 10 millions rays are 15 launched per model. For each ray, we measure the average, median 16 and maximal number of steps required to process the ray. Note that 17 the maximal number of rays are more important for SIMD architec-18 tures. Statitics are provided in Table 2 and 3. It is important to recall that our method relies on gradient computation for the evaluation 20 of directional derivative. As described in [ZBQC13], both field and 21 gradient can be computed for less than twice the evaluation cost of 22 a single field value. Furthermore, as both values are computed at 23 once, segment parameters are only fetched once from memory.

Figure 19 :

 19 Figure19: Individual times for each substep of the methods (occluding depth buffer, dexel buffer creation with and without occluders and ray processing). For models presenting a large number of depth layers, occluders provide non negligible improvement of rendering time.

Figure 20 :

 20 Fig.1

2

  new sampling strategy to existing rendering techniques for more 3 general implicit surfaces such as [GGPP20; Kee20] could provide 4 interesting optimizations.

Figure 21 :

 21 Figure 21: Comparison of grazing rays in a slice of the scene of Figure 11 (far left) for sphere-tracing, segment-tracing and our method. Note that homothetic distance primitives (Equation (7)) are used for this figure in order to allow direct comparison of the three methods.

  Processing interval (t begin , h b , t end , he). // h b , he are tuples containing Hermite data n ← 0, saved ← None while + + n < 32 do troot ← Oracle(t begin ,t end , i) if troot == ∞ and !saved then return f alse, ∞ else if troot < ∞ then hr ← HermiteData(troot ) if |hr.val| < ε then return true,troot end if b prim ← hr.prim > 0 and h b .prim < 0 if !saved and hr.val ≥ 0 and he.val < 0 and b prim then saved ← [(troot , hr), (t end , he)] end if if hr.val < 0 or b prim then t end , he ← troot , hr else t begin , h b ← troot , hr end if else if saved then (t begin , h b ), (t end , he) ← saved

	saved ← None
	end if
	end while
	return f alse, ∞
	Figure 9: A polynomial interpolation can present smaller varia-
	tions than the actual field function, resulting in possible missed iso-
	values. Until a root is isolated, we rely on rational Bézier curves to
	address this problem.

submitted to COMPUTER GRAPHICS Forum (2/2021).

Algorithm 1

Table 1 :

 1 Rendering Times (in milliseconds) averaged over full rotations around the object, at 2048x2048.

				UHD Graphics 630		nVidia Quadro P1000	nVidia GeForce 2080 RTX
				Ours	Sphere Sherstyuk	Ours	Sphere Sherstyuk	Ours Sphere Sherstyuk
		Fig.11(4)		290.4	551.9	565.9		84.2	300.4	262.5	7.9	28.2	18.4
		Fig.11(3)		284.4	374.4	526.2		80.0	169.0	230.3	7.9	16.1	17.0
		Fig.11(2)		243.0	388.0	384.1		127.8	302.0	268.8	12.0	27.2	22.5
		Fig.11(1)		169.5	244.6	289.6		74.6	169.0	182.1	7.3	15.6	15.3
		Fig.13(Middle)	121.3	204.9	178.4		31.8	120.0	79.9	6.2	12.9	7.7
		Fig.2		116.5	169.8	177.1		23.8	81.9	70.8	6.3	12.9	7.7
		Fig.15(Middle)	159.3	272.0	265.4		39.9	151.0	115.4	6.6	15.1	10.0
		Fig.1		129.7	247.8	202.4		33.6	152.7	93.3	6.7	22.0	13.2
		Fig.15(Left)	207.7	366.2	383.0		49.2	182.0	151.1	6.9	18.9	11.3
		Fig.15(Right)	117.5	187.1	181.1		29.1	106.8	76.3	5.9	10.4	7.0
		Fig.16		90.5	120.6	150.7		22.5	69.1	76.0	5.5	8.4	7.5
		Fig.12(Left)	294.5	390.0	447.8		119.4	236.0	246.0	11.3	25.6	22.7
		Fig.12(Middle)	-	-	-		-	-	-	62.5	259.9	97.3
		Intel UHD Graphics 630		nVidia Quadro P1000	nVidia GeForce 2080 RTX
	Average processing time (ms)	0 500 1,000			0 200 400				0 10 20		Occluders Dexel buffer DB (no occluders) Ours Sphere-tracing Sherstyuk's RT
		F i g .1	F i g .1 6 F i g .1 1 ( 2 ) F i g .1 2 ( L e f t )	F i g .1	F i g .1 6 F i g .1 1 ( 2 ) F i g .1 2 ( L e f t )	F i g .1	F i g .1 6 F i g .1 1 ( 2 ) F i g .1 2 ( L e f t )

Table 2 :

 2 submitted to COMPUTER GRAPHICS Forum (2/2021). Statistics calculated for our method and sphere tracing for Compact Polynomial. Similar results are calculated for the Cauchy kernel with decrease in median number of steps in the range between 36.3% and 97.9%

	Object	# of Skeletons	Our Method Avg Median Max Avg Median Sphere Tracing Max
	Fig.12(Left)	4722	5	4	37	16	11	24580
	Fig.11(1)	255	5	4	80	18	13	2389
	Fig.11(3)	201	6	6	84	19	14	2124
	Fig.15(Right)	127	5	5	27	36	24	2980
	Fig.1	213	5	5	132	68	22	42214
	Fig.15(Left)	55	5	5	23	73	48	8298
	Fig.11(4)	201	8	7	84	413	330	21770
	Fig.13(Middle)	95	5	4	23	31	23	4866
	Fig.15(Middle)	55	5	5	21	37	20	7559
	Fig.2	36	4	4	34	28	18	4372
	Fig.16	608	5	4	46	16	10	8345
	Object	# of Skeletons	Our Method Avg Median Max Avg Median Max Segment Tracing
	Fig.12(Left)	4722	5	4	53	9	7	383
	Fig.11(1)	255	5	5	55	9	8	549
	Fig.11(3)	201	5	4	24	9	7	340
	Fig.15(Right)	127	4	5	20	16	14	413
	Fig.1	213	5	4	65	15	11	736
	Fig.15(Left)	55	5	5	16	23	23	302
	Fig.11(4)	201	6	4	25	177	197	1680
	Fig.13(Middle)	95	4	4	18	15	14	277
	Fig.15(Middle)	55	4	6	18	13	16	531
	Fig.2	36	4	5	16	12	12	524
	Fig.16	608	5	6	74	8	7	348

Table 3 :

 3 Statistics calculated for segment tracing

  The Visual Computer 14 12.10 (Dec. 1996), 527-545. ISSN: 1432-2315. DOI: 10 . 1007 / 15 s003710050084. URL: https : / / doi . org / 10 . 1007 / Modeling System for Base Meshes of 3D Articulated Shapes". Com-19 puter Graphics Forum 29.7 (2010), 2169-2177. DOI: 10 . 1111 / j .

	16
		s003710050084 4.
	20
	21	1467-8659.2010.01805.x 1.
	22	[JTZ09] JIN, XIAOGANG, TAI, CHIEW-LAN, and ZHANG, HAILIN. "Im-
	23	plicit modeling from polygon soup using convolution". Vis. Com-
	24	put. 25.3 (Feb. 2009), 279-288. ISSN: 0178-2789. DOI: 10 . 1007 /
	25	s00371-008-0267-3 4.
	26	[KB89] KALRA, D. and BARR, A. H. "Guaranteed Ray Intersections with
		Implicit Surfaces". SIGGRAPH Comput. Graph. 23.3 (July 1989), 297-

2 Computer Graphics Forum 34.6 (May 2015), 198-210. DOI: 3 10 . 1111 / cgf . 12530. URL: https : / / hal . archives -4 ouvertes.fr/hal-01257198 5. 5 [GGPP20] GALIN, ERIC, GUÉRIN, ERIC, PARIS, AXEL, and PEYTAVIE, 6 ADRIEN. "Segment Tracing Using Local Lipschitz Bounds". Com-7 puter Graphics Forum (2020). URL: https : / / hal . archives -8 ouvertes.fr/hal-02507361 5, 13, 16. 9 [GPP*10] GOURMEL, OLIVIER, PAJOT, ANTHONY, PAULIN, MATHIAS, 10 et al. "Fitted BVH for Fast Raytracing of Metaballs". Computer Graph-11 ics Forum 29.2 (May 2010), xxx-xxx 4, 5. 12 [Har96] HART, JOHN C. "Sphere tracing: a geometric method for the 13 antialiased ray tracing of implicit surfaces".

17

[JLW10] JI, ZHONGPING, LIU, LIGANG, and WANG, YIGANG. "B-Mesh: 18 A

submitted to COMPUTER GRAPHICS Forum (2/2021).

ure 12) and random skeletons (see Figure11) using both low-end

Acknowledgements 6

This work was supported by the ANR IMPRIMA(ANR-18-CE46-

Such limitations have been tackled recently in the segment trac-16 ing algorithm [GGPP20], drastically decreasing the number of 0004). We thank Sylvain Lefebvre and Salim Perchy for their 8 proofreading. We thank Pierre-Alexandre Hugeron and Jonas Mar- We present here the worst case skeleton configuration used to compute the scaling factor of the sphere-cones' radii that is used for computing the occluding depth buffer.

Usage of the end-point correctors guarantees the existence of a minimal length of skeleton around any skeleton point q with radius τ (i.e. the length of the correctors on both sides of q) given that q is not part of an end-point corrector.

We want to find the minimal distance between the isosurface and the point q. Hence, we should study the field variation along arbitrary directions starting from q. For an arbitrary direction v, the minimal distance will be reached for the fastest decreasing field. Such field is obtained if the guaranteed length of skeleton is leaving q in the direction -v (kernels are decreasing functions). Similarly the influence of a given skeleton point is decreasing more rapidly for smaller prescribed radii.

Hence, the worst case scenario is reached for two segment primi-60 tives in the direction -v with radius τ in q and the fastest authorized 61 radius variation (see Figure inset). Thanks to the scale invariance of 62 the surface, we only need to compute the scaling factor for a unit 63 radius which we perform numerically (which corresponds to the 64 distance between q and the isosurface in the direction v -see Figure 65 inset). Note that this worst case configuration is actually indepen-66 dent from the actual skeleton configuration, no additional skeleton 67 processing is required.

Appendix B: Lipschitz constant computation 69

When applying sphere-tracing, computing the Lipschitz constant per primitive (e.g. ignoring the ray/primitive configuration) can have a large impact in presence of varying radius along primitives. Indeed, for the studied implicit surfaces the Lipschitz constant is inversely proportional to the smallest radius. In order to alleviate this problem in our sphere-tracing implementation, we compute the Lipschitz constant both per primitive and per ray. This can be done by computing the primitive point with the smallest radius whose kernel support is intersected by the ray. Based on Equation ( 18), the kernel support associated to a skeleton point is intersected by the ray if and only if:

where σ is the kernel scale. The smallest radius can then be deduced from these polynomial roots.


which separate from the residual input pulse due to group velocity walk-off. Also, because continuum contains such a broad spectrum that dispersion in the spectral wings becomes important, further propagation actually broadens the continuum significantly in time, rather than frequency. All the latter features, which only occur for longer (>1 cm) fiber lengths, are generally undesirable.

Previous measurements of the continuum had been limited to simple multishot spectral measurements, which have always yielded a broad, smooth and stable spectrum. In contrast, we performed much more powerful cross-correlation frequency-resolved-optical-gating (XFROG) measurements with a newly developed angle-dithered-crystal technique to achieve the huge bandwidth necessary to measure continuum. These measurements generate a spectrogram of the continuum, which visually displays its time and frequency characteristics, as well as yielding the first complete intensity and phase measurement of the continuum. Our measurements revealed the continuum intensity and color vs. time and yielded some very surprising conclusions. First, in strong contrast to previous simple multi-shot spectrometer measurements, they showed a very complex and unstable spectrum, which was in good agreement both with our theory and also with single-shot spectral measurements we made to confirm these unintuitive results. Thus, we found that the contin- O ne of the most exciting recent developments in optics has been the generation of ultrabroadband supercontinuum, accomplished simply by injecting readily available low-power ultrashort pulses into microstructure or tapered fibers. [2][3][4] These new fibers are nearly dispersion free at the Ti:sapphire laser wavelength, so the pulse remains short for a much longer distance than in conventional fibers (for many centimeters compared with a few hundred micrometers), dramatically increasing nonlinear optical effects. The resulting continuum's spectrum encompasses the entire visible and much of the infrared ranges, and spatially it is also highly coherent. Many applications, in areas ranging from metrology to medical imaging, have been proposed and demonstrated for this exotic light, but detailed measurement and understanding of this most complex pulse ever generated have both eluded researchers. Indeed, the incredible complexity of the supercontinuum provides a particular challenge to both models and measurement techniques.

In this work, we combined powerful new modeling capabilities and new measurement techniques to accurately model and measure the continuum's complete intensity and phase vs. time, revealing several complete surprises and important new strategies for using this fiber in the future.

In our simulations, for example, we found, quite surprisingly, that the extreme spectral width is actually achieved in the first centimeter of fiber. Also surprisingly, the strong spectral broadening in the first centimeter is also accompanied by strong pulse temporal compression.With further propagation, the already broad spectrum develops only asymmetry and complex temporal features, such as temporal pulse break-up, oscillations and distinct soliton pulses, uum spectrum is broad, but neither smooth nor stable! Also, they showed that the continuum pulse generated in long (many centimeter) fibers is quite longseveral picoseconds-in agreement with our simulations.

The figure shows the theoretical and experimental continuum spectrograms, revealing the complex structure and the typical spectral oscillations that occur with a long (16 cm) fiber. The oscillation frequency varies across the profile, with a mean frequency of about 115 THz.

These simulations and measurements clearly showed that, while the input pulse can propagate large distances in these fibers without distortion, the continuum cannot.Thus, the optimal approach to supercontinuum generation is to use a short, ~ 1 cm, fiber. Indeed, using such a fiber, we have recently succeeded in generating a supercontinuum pulse only 25 fs long-considerably shorter than the 40-fs pulse that created it-and also much smoother and much more stable. This short-fiber continuum is not only a nearly ideal pulse for most broadband applications, it is also potentially compressible to a few femtoseconds.

1 Rick

 1 Trebino, John Dudley and Xun Gu

Figure 1 .

 1 Figure 1. (a) Simulated supercontinuum spectrogram which shows the pulse break up and complex spectrum after propagating 16 cm. (b) M easured and retrieved (a check on the measurement) continuum spectrograms showing the complex behavior predicted. In both cases, the intensity is shown below and the spectrum at right. (The original Optics Express paper contained a movie that showed the continuum spectrogram vs. distance along the fiber.)


cleosomes, which is valid as soon as the internucleosomal distances are larger than the interaction range; in any cases, the model describes the geometric contribution to chromatin elastic behavior originating from the DNA elastic properties.

Our model, quite similar to the so-called two-angle model of Woodcock et al. [START_REF] Woodcock | A chromatin coiling model that incorporates linker variability generates fibers resembling the native structures[END_REF], amounts to build the chromatin fiber nucleosome after nucleosome. The rotational positioning of a nucleosome with respect to the previous one is entirely prescribed by the twist angle of the DNA linker connecting them. In the relaxed state, it is equivalently prescribed by the linker length n (in bp) since the twist angle is then precisely equal to τ = 2π n/10.6 (recall that there is 10.6 bp per DNA turn). Each step of the assembly is described analytically, but we implemented this model within a Maple program in order to handle an arbitrary number of nucleosomes and to perform a quantitative analysis of the chromatin fiber geometric properties.

Structural properties

The implementation of the above model straightforwardly gives all the geometric properties of the chromatin fiber. Homogeneity of the local architecture (uniform values for n and Φ) enforces a helical symmetric structure for the fiber; we shall speak of "chromatin superhelix" (SH), and denote ⃗ A its axis. A complete study of the symmetry properties and quantities relevant to the mechanics of the fiber has been presented in a previous paper [START_REF] Ben-Haïm | Chromatin: a tunable spring at work inside chromosomes[END_REF]. We here only underline the high structural sensitivity of the fiber, evidenced on Figure 1: for various values of n and Φ, we obtain the whole variety of structures suggested from experimental studies [1] [12] [13] [START_REF] Widom | Structure, dynamics and function of chromatin in vitro[END_REF]. These qualitative results support our minimal two-parameter modeling, since it is sufficient to recover the actual diversity of chromatin structures. A noticeable point, evidenced on Figure 2, is the fact that the diameter remains about 30 nm whatever the structure; the so-called 30nm-fiber observed experimentally might thus correspond to far different structures. Excluded volume has to be taken into account to discard unrealistic structures; we determined a sufficient condition to avoid steric hindrance, which is shown on Figure 2.

One of the conclusions of our structural study is the fact that the connection between the microscopic parameters and the SH geometric characteristics is too complex and multivariate to get small-scale informations from structural observations of the fiber; we now turn to the analysis of its mechanical properties. Figure 2: Diameter of the chromatin fiber (upper blue curve) and its pitch P (lower red curve). The thin green line gives a lower bound Pc for the pitch above which no steric hidrance occurs within the fiber, in constrast to the dashed regions where P < Pc.

Elastic properties

In order to obtain the leading order of the elastic constants of the fiber, we merely consider a straight fiber with cylindrical symmetry. Describing the elastic properties of the chromatin fiber is thus nothing but a problem of spring mechanics [START_REF] Love | A treatise on the mathematical theory of elasticity, 4th edition, chapters XVIII and XIX[END_REF]. Nevertheless, the architecture of this "spring" is much more complex than a simple helical coiling and we expect that the detailed structural features of the chromatin assembly strongly influence the behavior at the fiber scale. Our modeling allows us to compute analytically the elastic constants describing the linear response of the fiber when a force ⃗ F and a torque ⃗ M are applied at its ends, hence to investigate quantitatively how the elastic response of the fiber varies with its relaxed structure. Note that this structure is itself controlled by the "microscopic" structural parameters n and Φ.

Our computation lays on the identification of the elastic energy of the fiber when a force ⃗ F and a torque ⃗ M are applied at its ends and the sum of the elastic energies of its linkers (bend and twist energies) in such a constrained configuration. This identity follows from the absence of internucleosomal interactions and to our focus on the linear response regime (no modification of the fiber architecture at this order); it amounts to restrict to the elastic contribution originating from linker elasticity within the relaxed chromatin structure. The first step is then to relate the "macroscopic" stresses ( ⃗ F , ⃗ M ) to the distribution of stresses along linker DNA. We denote ⃗ fj (s) the force and ⃗ mj(s) the torque exerted at the point Pj(s) of linker j (with arclength s) by the upstream part of the fiber. We focus on the universal, rotationally symmetric behaviour of the fiber. Writing standard equilibrium equations of spring mechanics, and using the linear response hypothesis to simply sum up the effects of the different stresses applied to the fiber (stretching force, bending torque and torsional torque), we obtain:

⃗ fj (s) = ⃗ F ⃗ mj (s) = ⃗ M -[ ⃗ Oj (s)Pj(s)] ∧ ⃗ F (1)
where Oj (s) denotes the orthogonal projection of Pj (s) onto the SH axis ⃗ A. It gives the relation between the global stresses ⃗ F and ⃗ M exerted at the fiber ends and the local stresses ⃗ fj (s) and ⃗ mj(s) experienced at the linker level, at each point of the DNA path. The term [ ⃗ Oj (s)Pj(s)] ∧ ⃗ F reflects the involvement of the fiber architecture in the expression of the local torque ⃗ mj(s). This set of equations is the core of our analysis since it allows to perform explicitly the change of description level, namely passing from the DNA level to the fiber level and reciprocally.

As we restrict to the description of harmonic elasticity, the coefficients of ⃗ F and ⃗ M in the elastic energies of the linker will be computed within the relaxed SH. We underline that we do not need to compute the constrained shape of the linkers to describe the linear response of the fiber to applied force and torque. The cylindrical symmetry of the fiber ensures that ⃗ F is directed along the SH axis ⃗ A at equilibrium. ⃗ M decomposes into a bending torque ⃗ M b and a torsional torque

Mt according to ⃗ M = ⃗ M b + Mt ⃗ A.
The well-established elastic energy of a linker [3] [4] can then be expressed as a function of F , Mt and M b thanks to (1); it yields a quadratic expression

E linker = asF 2 + atM 2 t + 2atsMtF + a b M 2 b ( 2
)
which is to be identified with the elastic energy E fiber (F, Mt, M b ) of a length D of chromatin fiber, where D is the distance between two successive nucleosomes measured along the SH axis (projected distance).

Plugging in the quadratic expression E linker (F, Mt, M b ) the linear response ansatz relating the stresses F , Mt and M b experienced by the fiber and its strains u (relative extension), Ω (twist rate) and ρ (local curvature):

F Mt M b = γ kBT g 0 kBT g kBT C 0 0 0 kBT A u Ω ρ (3)
naturally embeds the fiber and its elastic behaviour within a current continuous model known as the "extensible worm-like rope" ( EWLR) [START_REF] Marko | DNA under tension: overstretching, undertwisting and relaxation dynamics[END_REF], in which the density of elastic energy writes, as a function of the canonical variables u, Ω and ρ:

ϵEWLR = kBT Aϱ 2 2 + kBT CΩ 2 2 + γ u 2 2 + kBT g Ωu (4)
A is the bend persistence length of the fiber, C its twist persistence length, γ its stretch modulus (dimension of a force) and g the twist-stretch coupling constant (no dimension). Note the expected vanishing of the twist-bend and stretch-bend coupling (for symmetry reasons [START_REF] Marko | DNA under tension: overstretching, undertwisting and relaxation dynamics[END_REF]. We may now express the elastic constants of the fiber in terms of the coefficients a b , at, ats and as, obtained as functions of the DNA elastic constants and fiber geometry, by a term-wise identification of the coefficients of u 2 , Ω 2 , uΩ and ρ 2 in the following identity:

E fiber ≡ D kBT Aϱ 2 2 + kBT CΩ 2 2 + γ u 2 2 + kBT g Ωu = E linker [F (u, Ω), Mt(u, Ω), M b (ρ)] (5)
from which follows analytical expression [START_REF] Ben-Haïm | Chromatin: a tunable spring at work inside chromosomes[END_REF] of the elastic constants A, C, g and γ. We present on Figure 3 the results for Φ = 50 o (high salt concentration, with histone H1, i.e. physiological conditions)

when the linker length is varied. The curves evidence a strong sensitivity of elastic coefficients with respect to the fiber structure, here controlled by n. This sensitivity originates in the sensitivity of the SH pitch: accumulation of turns in a spring leads to a highly flexible accordion-like behavior. It is to note that the persistence length A and C are almost identical, both close to the SH pitch and smaller than the corresponding persistence lengths A and C of DNA. More strinkingly, the stretch modulus of the chromatin fiber, around 5 pN, is to be compared to that of DNA (γDNA ≈ 1100 pN): whereas DNA can be seen as non extensible, chromatin is highly extensible; in particular, entropic and enthalpic elastic regimes (well-separated in the case of DNA) will overlap. We finally underline a key feature: the twist-stretch coupling constant, reflecting the chirality of the chromatin fiber, exhibits steep variations and sign changes when the linker length is varied. The same twisting of the fiber will either condense or decondense it according to the relaxed structure.

Experimental implications

A first issue involving fiber elasticity is to describe the linear elastic response of the fiber to global stresses, i.e. a force and a torque applied at its ends. This issue refers to micromanipulations in which a single chromatin fiber is pulled [START_REF] Cui | Pulling a single chromatin fiber reveals the forces that maintain its higher-order structure[END_REF] (and possibly will be twisted). Describing the fiber within the EWLR model, yet thoroughly investigated in the context of DNA by Marko [START_REF] Marko | DNA under tension: overstretching, undertwisting and relaxation dynamics[END_REF], allows an interpretation of the force-extension curves in terms of geometric and mechanical parameters of the fiber. Experimental results of Cui and Bustamante [START_REF] Cui | Pulling a single chromatin fiber reveals the forces that maintain its higher-order structure[END_REF] are in good agreement with our predictions (note that their measure rests on relaxation curves, for which there is no interaction between nucleosomes, consistently with our modeling). Indeed, they measured a bend persistence length A ≈ 30 nm and a stretch modulus γ ≈ 5 pN, which are obtained simultaneously in our modeling for n ≈ 42, or any integer multiple of DNA pitch; we there recover the nucleosome phasing observed in biological experiments [START_REF] Yao | Twist constraints on linker DNA in the 30-nm chromatin fiber: implications for nucleosome phasing[END_REF].

Our study should now be supplemented by investigating the effect of a dispersion in the values of n and Φ. It is likely that the dispersion will induce a kind of averaging. Continuity arguments lead us to expect the sensitivity of the structural and elastic properties to remain -although possibly smoothed out -together with its possible relevance in biological regulation. Work is in progress on this question.

To circumvent this limitation, we suggest to take advantage of the interplay between linker intercalation and the chromatin fiber mechanics. The twist change ∆τ achieved in the linker by intercalator molecules (for instance with ethidium bromide in vitro) amounts to shift the linker length by ∆n = 10.6 × ∆τ /2π. Plotting the various elastic constants as a function of ∆n should provide similar curves as those shown on Figure 3, however smoothed out. We therefore suggest that it should be possible to scan experimentally these curves by varying continuously the concentration of intercalators added in the buffer when performing pulling fiber experiments. (a) Bend persistence length A (solid red line) in nm and stretch modulus in pN (dashed blue line); here is shown the effective modulus γ ef f = γ-kBT g 2 /C which is involved in force-extension curves [START_REF] Marko | DNA under tension: overstretching, undertwisting and relaxation dynamics[END_REF] [11] hence directly measurable. Our computations predict that A and the twist persistence length C are almost equal, both close to the pitch P of the SH (Figure 2). Experimental measurements of Cui and Bustamante [START_REF] Cui | Pulling a single chromatin fiber reveals the forces that maintain its higher-order structure[END_REF] yield a value about 30 nm for A and 5 pN for γ ef f (horizontal lines).

(b) Twist-stretch coupling constant g (no dimension). Note that the vanishing of g corresponds to a change of chirality for the fiber.

Biological relevance

Another issue, directly relevant to the in vivo functioning of chromatin, is to investigate both the effect of DNA-protein binding on the fiber structure, and the effect of fiber stresses on the linker DNA-protein binding. Our theoretical approach gives a framework to such studies by relating the structure and mechanics at the DNA scale and those at fiber scale.

On one hand, our approach allows to describe the response of the fiber to local, internal stresses as those created by intercalators, groove-binding proteins, or any induced change in the fiber assembly at the linker scale. Mechanical sensitivity is likely to provide efficient switches for processes involving a conformational transition of the fiber, for instance the decondensation required prior transcription or, at the opposite, fiber compaction currently associated with gene silencing [START_REF] Widom | Chromosome structure and gene regulation[END_REF]. A possible tuning mechanism lays on the variation of Φ, controlled in particular by the presence of linker histone, salt concentration and histone tails binding affinities.

On the other hand, protein-DNA interactions involved in gene-regulated subcellular events take place within chromatin fiber. In consequence, they are likely to be influenced, if not directly regulated, by chromatin structure and by the stresses generated in the linkers when the chromatin fiber is somehow constrained. In particular, we expect such an involvement of chromatin structure and multilevel mechanical properties to be at work in epigenetics and gene imprinting [19].

Conclusion

The chromatin scale is precisely the scale of nanomechanics: at this scale, we expect a strong and direct interplay between the biological functioning, monitored by various enzymes, and the mechanical properties of the chromatin fiber. Our complete and quantitative analysis gives clues to discuss how the elastic properties of the chromatin fiber might at the same time favour DNA compaction into the chromosomes and allow local and controlled decondensation of chromatin involved in gene expression.

Figure 1 :

 1 Figure 1: Varying the linker length n and the entry/exit angle Φ recovers all the structural diversity observed experimentally, for instance ribbon (Φ = 90 o , n = 58 bp), columnar structure (Φ = 50 o , n = 42 bp) or double-helical structure (Φ = 60 o , n = 56 bp).

Figure 3 :

 3 Figure 3: Elastic constants of the chromatin fiber versus linker length n (in bp) for Φ = 50 o . A similar sensitivity when n varies is observed for other values of Φ, or when varying Φ at fixed n.(a) Bend persistence length A (solid red line) in nm and stretch modulus in pN (dashed blue line); here is shown the effective modulus γ ef f = γ-kBT g 2 /C which is involved in force-extension curves[START_REF] Marko | DNA under tension: overstretching, undertwisting and relaxation dynamics[END_REF] [11] hence directly measurable. Our computations predict that A and the twist persistence length C are almost equal, both close to the pitch P of the SH (Figure2). Experimental measurements of Cui and Bustamante[START_REF] Cui | Pulling a single chromatin fiber reveals the forces that maintain its higher-order structure[END_REF] yield a value about 30 nm for A and 5 pN for γ ef f (horizontal lines).