extfieldsofstudy:
Computer Science
Mathematics

s2fieldsofstudy:
Computer Science
Mathematics

text:
Coded Merkle Tree: Solving Data Availability Attacks in Blockchains

In this paper, we propose coded Merkle tree (CMT), a novel hash accumulator that offers a constant-cost protection against data availability attacks in blockchains, even if the majority of the network nodes are malicious. A CMT is constructed using a family of sparse erasure codes on each layer, and is recovered by iteratively applying a peeling-decoding technique that enables a compact proof for data availability attack on any layer. Our algorithm enables any node to verify the full availability of any data block generated by the system by just downloading a $\Theta(1)$ byte block hash commitment and randomly sampling $\Theta(\log b)$ bytes, where $b$ is the size of the data block. With the help of only one connected honest node in the system, our method also allows any node to verify any tampering of the coded Merkle tree by just downloading $\Theta(\log b)$ bytes. We provide a modular library for CMT in Rust and Python and demonstrate its efficacy inside the Parity Bitcoin client.

Introduction
Blockchains (e.g., Bitcoin [25] and Ethereum [34]) maintain a ledger of ordered transactions, organized into a chain of blocks. Starting from the genesis block, network nodes extend the ledger by creating and appending more blocks, following specific block generation rules (e.g., the longest-chain rule is used in Bitcoin [25]). The transactions in the received blocks are validated by full nodes which download the entire block tree. However, for better scalability, it is imperative for a blockchain to allow light nodes, which may only be interested in verifying some specific transactions.
In Bitcoin [25,1], light nodes are implemented using the Simple Payment Verification (SPV) technique: a Merkle tree is constructed for each block using the transactions as the leaf nodes, and the Merkle root is stored in the block header. Utilizing the Merkle root, a light node can verify the inclusion of any transaction in a block through a Merkle proof. Light nodes and SPV have been leveraged extensively to scale computation and storage of blockchain systems over resource-limited nodes (e.g., smartphones) [3,7,4,5,24,35,17,15,18].
Besides inclusion, what is more important for a light node is to validate the transaction based on the ledger state. Due to limited resources, a light node cannot download the entire ledger. Instead, it could use the depth of the block that contains this transaction as a proxy. That is, the deeper this block is buried into the chain, the more confident the light node is about the validity of the transaction. However, for it to work, a majority (in terms of hashing power, stakes, etc.) of full nodes must be honest and must follow protocol. Further, there is a significant tradeoff between confirmation latency (due to the depth) and the security about transaction validity.
Therefore, efforts to study 1) the scenario of a light node being connected to dishonest majority of full nodes, and 2) how to achieve faster confirmation at light nodes are becoming a major research direction [25,9,16,10]. The overall idea is to design new block structures that allow full nodes to generate and broadcast succinct fraud proofs of individual transactions. This way, a light node will be able to timely verify fraud transactions and blocks as long as it is connect to at least one honest full node. One efficient construction that utilizes the roots of the intermediate state Merkle trees after executing a subset of transactions is proposed in [9]. However, it is vulnerable to the so-called "data availability attack" described in [9], for which [9] proposed an erasure code based solution. Stating the data availability attack formally and solving it comprehensively is the main goal of this paper. Data availability attack. A malicious block producer 1) publishes a block header, so that light nodes can check transaction inclusion; but 2) withholds a portion of the block (e.g., invalid transactions), so that it is impossible for honest full nodes to validate the block and generate the fraud proof.
Although the honest full nodes are aware of the data unavailability, there is no good way to prove it. The best they can do is to raise an alarm without a proof. However, this is problematic because the malicious block producer can release the hidden parts after hearing the alarm. Due to network latency, other nodes may receive the missing part before receiving the alarm and, thus, cannot distinguish who is prevaricating. Due to this, there is no reward and punishment mechanism that can properly reward honest full nodes while also deterring false alarms and denial-of-service attacks.
Therefore, for fraud proofs to work, light nodes must determine data availability by themselves. This leads to the following key question: when a light node receives the header of some block, how can it verify that the content of that block is available to the network by downloading the least possible portion of the block? Need to encode the block. Since a transaction is much smaller than a block, a malicious block producer only needs to hide a very small portion of a block. Such hiding can hardly be detected by light nodes unless the entire block is downloaded. However, by adding redundancy to the data through appropriate erasure codes [21], any small hiding on the origin block will be equivalent to making a significant portion of the coded block unavailable, which can be detected by light nodes through randomly sampling the coded block with exponentially increasing probability. As a counter measure, a malicious block producer could instead conduct coding incorrectly to prevent correct decoding. Light nodes rely on honest full nodes to detect such attacks and prove it through an incorrect-coding proof.
For example, an (n, k) Reed-Solomon (1D-RS) code [30] encodes k data symbols into n coded symbols, and any k out of these n coded symbols can be used to decode the k data symbols. Thus, to prevent decoding, a malicious block producer will have to make at least n−k+1 coded symbols unavailable, which yields a detection probability of 1 − (1 − k/n) s after sampling s distinct coded symbols uniformly at random. But an incorrect-coding proof will consist of k coded symbols, which is of the same size as the original block and thus is too large. This cost is alleviated to √ k in [9] by using two-dimensional RS codes (2D-RS), at the costs of reduced sampling efficiency, and increased block hash commitments of 2 √ n Merkle roots to verify the coding correctness within each dimension. In addition, 1D-RS and 2D-RS codes have a high decoding complexity of O(k 2 ) and O(k 1.5 ), respectively.
In summary, with erasure coding, a light node pays 3 download costs for data availability, including block hash commitments, symbol sampling, and incorrectcoding proofs. Among them, the incorrect-coding proof cost must be minimized to defend fake proofs, for which both 1D-RS and 2D-RS are sub-optimal. Our contributions. In this paper, we propose SPAR (SParse frAud pRotection), the first data availability solution that promises order-optimal performance on all the metrics, including 1) a constant block hash commitment size; 2) a constant sampling cost for a given confidence level on data availability; 3) a constant incorrectly-coding proof size; and 4) linear decoding complexity (Table 1).
At the core of SPAR is a novel cryptographic hash accumulator called coded Merkle tree (CMT). Starting from the bottom, CMT iteratively encodes layers of the tree and uses the hashes of the coded layer as the data for the next layer. A light node can detect the availability of the entire tree through the Merkle proofs of bottom layer leaves. With the entire tree available, SPAR uses a novel hashaware peeling decoder and a special ensemble of random LDPC (low-density parity-check) codes to maximize sampling efficiency, minimize incorrect-coding proof to one parity equation, and achieves linear decoding complexity. SPAR and CMT implementation. We have developed a complete and modular CMT library in Rust and Python [2]. We have also implemented SPAR in the Bitcoin Parity client [33], which outperforms state of the art [9] by more than 10-fold in hash commitments, incorrect coding proof, and decoding speed. Fig. 1: The connections between full nodes and light nodes. Every honest (green) node is connected to at least one honest full (square) node. A light (circle) node may connect to more malicious (red) full nodes than honest full nodes.
Related works. This work was inspired by pioneering research in [9], which proposes succinct fraud proofs and 2D-RS based data availability solution. Besides this work, coding also improves scalability of blockchains in other areas: [27] studies the coding efficiency of distributed storage systems [8,14,13,29]. In a related vein, [28] uses a combination of Shamir's secret sharing [31] (for storing the headers and the encryption keys) and private distributed storage (for the blocks) to reduce the storage overhead while guaranteeing data integrity and confidentiality. [20] uses Lagrange coding to simultaneously scale storage, computation, and security in a sharded blockchain [23,19], via cross-shard coding.

Security Model
The core functionality of a blockchain is to produce, verify, and accept/store valid data blocks in a consistent but decentralized manner. A data block, denoted by B, is a byte string of length b that carries a batch of transactions. B is valid for acceptance if and only if every single transaction in it is valid (e.g., enough balance, no double spending). Thus incomplete data blocks are tantamount to being unacceptable. Data incompleteness is not a threat to a node that fully downloads the block. However, state-of-the-art blockchain systems also run light nodes which do not download the blocks in entirety. We next describe these two types of nodes formally (see Fig. 1).
Full nodes are able to produce blocks (e.g., by batching submitted transactions), and to download and verify blocks produced by other full nodes. Upon acceptance, they store the entire block locally. Upon rejection, they broadcast a fraud proof to alert the other nodes. We note, however, that malicious full nodes do not necessarily follow such requirements, and can act arbitrarily.
Light nodes can only afford to download a small amount of data from each block and perform simple computations such as hash checks and fraud proof verification, but not to operate on whole data blocks. By accepting a block B, they only store its hash commitment D = g(B). Here g() is a hash accumulator such as Merkle tree generator, which will allow it to use D to verify the inclusion of any transaction in B through a Merkle proof. Without loss of generality, we assume light nodes are honest, as they are not able to deceive full nodes.
We assume the following network model: 1. Reliable communication: Two directly connected nodes can reliably communicate via both unicast and broadcast without message loss or corruption.
2. Connectivity: Every honest (full and light) node is directly connected to at least one honest full node. In other words, we assume a connected sub-graph of honest full nodes, and all the light nodes are connected to it.
3. The network is synchronous. These three assumptions together means that a valid message sent by an honest node can be received by all appropriate honest nodes (e.g. blocks for honest full nodes, block headers and fraud proofs for light node) within a fixed delay if every honest node re-broadcasts it.
4. The network allows nodes to send messages anonymously. Importantly, our network model allows dishonest majority, i.e., each light node can be directly connected to more malicious full nodes than honest ones. Due to this, a light node cannot determine the completeness of a block through its connected full nodes, via a majority vote for instance.
A malicious block producer is thus motivated to conduct a data availability attack, where it 1) does not fully disclose B, so that honest full nodes are not able to verify B; and 2) broadcasts D, so that itself and its colluding full nodes can forward D to their connected light nodes and deceive them that the B that satisfies g(B) = D is valid for accepting. Thus, the key for a light node to protect itself from accepting a fraudulent block is to make sure that B is fully available. This gives rise to the main problem we try to address in this paper: Data availability problem: Upon receiving a hash commitment D, how can a light node efficiently verify that a data block B that satisfies g(B) = D is fully available to the system? A simple strategy for a light node is to randomly sample portions of B, and determine that it is unavailable if it does not receive all requested portions. Since the size of a transaction is usually much smaller than the block, a malicious block producer only needs to hide a very small portion (say, e.g., a few hundred bytes) of a fraudulent block, which can hardly be detected through random sampling.
A malicious block producer could also conduct selective disclosure: when requested by light nodes, the malicious block producer may select a subset of the light nodes and fully disclose their requested portions, as long as the total disclosed portions do not reveal B. These light nodes will be deceived about the availability of B and will accept it, as no fraud proof of B can be produced. Thus, as similarly done in [9], we characterize the security of the above described system using the following measures: Soundness: If a light node has determined that a data block is fully available, then at least one honest full node will be able to fully recover this data block within a constant delay.
Agreement: If a light node has determined that a data block is fully available, then all the other light nodes in the system will determine that the data block is fully available within a constant delay.
Recently, an erasure coding-assisted approach was proposed in [9] to improve sampling efficiency and suppress the data availability attack. In the next section, we will motivate this approach and overview the challenges it faces.

Overview of Erasure Coding Assisted Approach
An (n, k) erasure code evenly partitions a block B of b bytes into k data symbols of b k bytes each as B = [m 1 , · · · , m k ], and linearly combines them to generate a coded block with n > k coded symbols, C = [c 1 , · · · , c n ]. The ratio r = k/n is called the coding rate. The n hashes of these coded symbols are accumulated to obtain the hash commitment D of C, which is published with C. With a good erasure code, a block producer's hiding of one data symbol is equivalent to making the value of many coded symbols unavailable to the system. In general, a pair of good erasure code and decoding algorithm yields a large undecodable ratio α, which is the minimum fraction of coded symbols a malicious block producer needs to make unavailable to prevent full decoding. Such hiding can be caught by a light node with an exponentially increasing probability of 1−(1−α) s through randomly sample s coded symbols when n is large, indicating that O(1) samples are sufficient. Below is an example. Example 1. Uncoded v.s. coded sampling efficiency. Given a block of 4 data symbols [m 0 , · · · , m 3 ], a block producer generates 8 coded symbols as follows: (1) To prevent decoding through hiding, a malicious block producer must either publish no more than 3 data symbols or no more than 5 coded symbols. Both will make at least 3 coded symbols unavailable to the system (α = 3 8 ). Such unavailability can be caught with a probability of 1 − 5 8 · 4 7 = 64.3% after randomly sampling 2 distinct coded symbols. In contrast, without coding, the hiding of one data symbol can be caught with a probability of only 1 − 3 4 · 2 3 = 50%.
To counter erasure coding assisted random sampling, a malicious block producer could conduct an incorrect-coding attack : It generates coded symbols that fail the parity equations (the equations describing the linear relations between coded symbols in Example 1) specified by the erasure code, and generates D using these invalid coded symbols. This way, it can pass light node random sampling through hiding only one data symbol and publishing most of the coded symbols, which will not allow honest full nodes to correctly decode B.
Fortunately, this attack can be detected by honest full nodes by comparing the decoded block with the commitment D. Upon detection, an honest full node can generate an incorrect-coding proof, which consists of the coded symbols of failed parity equation(s) and appropriate hash commitments, so that light nodes can verify them and reject the block. Using Example 1, an incorrect coding proof about c 4 = c 0 + c 1 could be c 0 and c 1 with matching Merkle proofs, plus the Merkle proof of c 4 , which, however, does not match the value of c 0 + c 1 .
To keep incorrect coding proofs small, [9] applies 2D-RS (2-dimensional Reed-Solomon) code. The k data symbols are placed as a √ k × √ k square, then a ( √ n, √ k) RS code is applied to every row/column. The resulted 2 √ n rows/columns yield 2 √ n Merkle roots, which are downloaded by light nodes as block header. Each root allows a light node to verify the associated row/column by decoding it using any √ k coded symbols of it (from incorrect-coding proof) and reproducing the root. Thus, 2D-RS offers light nodes 1) a header cost of Here log b is due to logarithmic growth of Merkle proof size with b. In this paper, we propose SPAR (SParse frAud pRotection), the first solution to the data-availability problem that is order-optimal in all the above three metrics: a header cost of O(1), a sampling cost of O(log b), and an incorrectcoding proof size of O(log b). To this end, SPAR leverages four core components: 1. a novel hash accumulator named coded Merkle tree (CMT), which encodes every layer of the tree to protect the availability of the entire tree. This way, the Merkle proof of every coded symbol will be available, which will enable every parity equation to be committed and verified alone; 2. a dedicated sampling mechanism that enables a light node to check the availability of the entire CMT by sampling O(log b) bytes plus one Merkle root; 3. a hash-aware decoding algorithm that is able to detect and prove any single failed parity equation, provided the Merkle proofs of all the coded symbols; 4. a special ensemble of random LDPC (low-density parity check) codes with a constant parity equation size and a constant undecodable ratio under the above hash-aware decoding algorithm, which protects all CMT layers equally.

Detailed Description of SPAR
In this section, we describe the four core components of SPAR: the construction of the coded Merkle tree by the (honest) block producer, the sampling mechanism of the light nodes, the decoding and alerting operations of the honest full nodes, and the erasure codes used by SPAR. At the end of this section, we will summarize the action space of each node in the network.

Construction of coded Merkle tree
In SPAR, an honest full node detects and proves incorrect-coding using the membership proofs of all the d coded symbols in one parity equation and the values of at least d − 1 of these coded symbols. Since any parity equation can be compromised, a light node needs to make sure the membership proofs of all the n coded symbols are available at honest full nodes. In other words, it needs to make sure the entire Merkle tree is available.
To this end, we propose CMT. At a high level, CMT applies erasure coding to every layer of the tree, where the data symbols of a layer are generated using the hashes of the coded symbols of its child layer. This way, a light node can check the availability of every layer through random samplings, whilst an honest full node can detect and prove the incorrect coding at any layer, with the help of the hashes of this layer provided at its parent layer. Fig. 2: Coded Merkle tree for k = 16 data symbols using rate r = 1 2 erasure codes. Each data symbol of a higher layer is constructed by batching the hashes of qr = 2 data symbols and q(1 − r) = 2 parity symbols of its child layer.
More specifically, given a block of k data symbols, a rate-r (r ≤ 1) systematic erasure code with an undecodable ratio of α is applied to generate n = k/r coded symbols, where the first k are the original data symbols and the remaining n − k are called parity symbols (hence the name systematic). Then the hashes of every q coded symbols are batched as one data symbol for the next (parent) layer. This yields a total of n/q data symbols for the next layer, which will be encoded using a smaller (in terms of k) rate-r systematic code with the same undecodable ratio. This iterative encoding and batching process stops once there are only t (t 1) hashes in a layer. These t hashes are the root of the CMT, and will be included in the block header and published with the original data block.
CMT layer size reduces at a rate of qr. Thus, qr > 1 for CMT to converge. In addition, to enable efficient sampling of both data and parity symbols (will discuss next), batching is interleaved, namely, the q coded symbols whose hashes are batched together consist of qr data symbols and q(1 − r) parity symbols. An example of CMT with k = 16, r = 1 2 , q = 4, and t = 4 is illustrated in Fig. 2. Indeed, a classic Merkle tree is a special CMT with r = 1 and q = 2.

Sampling Mechanism of Light Nodes
In SPAR, a light node randomly samples the base layer coded symbols with their Merkle proofs to decide the availability of the base layer. The special structure of CMT allows them to further utilize these proofs to efficiently sample higher layer symbols to decide the availability of higher layers.
Similar to a classic Merkle tree, the Merkle proof of a base layer symbol in CMT consists of all the sibling hashes between this symbol and the root. The only difference is that the number of sibling hashes per layer is now q − 1 instead of 1, which effectively provides the light node one data symbol from every intermediate layer. Thus, when a light node randomly samples s distinct base layer coded symbols, the associated Merkle proofs will automatically sample, at no extra cost, s distinct data symbols from every intermediate layer w.h.p.
To properly check the availability of an intermediate layer, a light node should also randomly sample about (1 − r)s parity symbols from this layer. To avoid downloading extra Merkle proofs for these parity symbols and to minimize the correlation between the samplings 5 , SPAR samples parity symbols of intermedi- For the 2 dash-circled yellow parity symbols, with a probability of 1 − r = 0.5, one of them (with equal chance) will be included in the proof for sampling purpose. So will the 2 dash-circled blue parity symbols. ate layers probabilistically: For every pair of parent and child intermediate layer, if a parent layer data symbol is sampled, then with probability 1 − r, one of its q(1 − r) child parity symbols (thanks to interleaved batching) will be sampled uniformly at random. Thus, the response size of one sampling request will be: where b k is the base layer symbol size, y is the hash size (e.g., 32 bytes), y(q −1) is the size of the partial data symbol from an intermediate layer for Merkle proof, yq(1 − r) is the average size of probabilistically sampled parity symbol from an intermediate layer, and log qr k rt is the number of layers. See Fig. 3 for sampling response of a coded symbol on the based layer of the CMT in Fig. 2.
Finally, to counter selective disclosure conducted by the malicious block producer, a light node will make the s requests separately, anonymously, with replacement, and with some delay between every two requests. This will prevent the malicious block producer from selectively deceiving any particular light node, or deceiving the set of light nodes that make requests at the beginning. Therefore, every light node will have the same chance to catch data availability attack.

Hash-aware peeling decoder and incorrect-coding proof
A hash-aware peeling decoder is similar to conventional LDPC peeling decoder. Given the hashes of all the n coded symbols and (1 − α)n coded symbols of a layer, it iteratively solves degree-1 parity equations and check each decoded symbol against its hash and associated parity equations (Algorithm 1). This way, the detection and proof of incorrect-coding is minimized to one parity equation.
The key condition for the peeling decoder to work is that the hashes of all the coded symbols are available. This is assured by CMT: By first downloading the root, the decoder will have all the hashes needed to decode the previous layer. Once this layer is successfully decoded, the decoded data symbols will provide all the hashes needed to decode its child layer. This top-down decoding continues Algorithm 1: Hash-aware peeling decoding algorithm 1 Inputs: the hashes of all n coded symbols and (1 − α)n coded symbols; 2 Initial check: checks all the degree-0 parity equations (i.e., those whose coded symbols are all known). If any parity equation is failed, report an incorrect-coding proof and exit; 3 while not all the k data symbols are recovered do 4 Find a degree-1 parity equation, which only has one unknown coded symbol; 5 Recover this coded symbol and verify it with its hash. If failed, report an incorrect-coding proof and exit; 6 Check all the associated degree-0 parity equations. If any parity equation is failed, report an incorrect-coding proof and exit. 7 end until the data block is decoded, or incorrect-coding is detected at one of the layers. To prove a failed parity equation that consists of d coded symbols, the decoder only needs to provide the Merkle proofs of these coded symbols, and the value of d − 1 coded symbols. Note that the higher the failed layer, the shorter the Merkle proof of each symbol in the incorrect-coding proof.
In addition, the peeling decoder only works if 1) there are (1 − α)n coded symbols available, and 2) that these coded symbols allow the recovery of all the k data symbols. While the first condition is checked by light nodes through random sampling, the second condition requires us to find, for every layer, a erasure code whose undecodable ratio is α under peeling decoding. The best performance is achieved if the codes are extremely sparse (with a small d) and have a large α. We now present such an ensemble of LDPC codes.

Construction of Erasure Code
An (n, k) erasure code can be described by an n × (n − k) parity check matrix H, where each column of H describes a parity equation, such that CH = − → 0 for any valid codeword C. In addition, every stopping set of H corresponds to a set of coded symbols whose hiding will prevent the full recovery of data symbols using peeling decoder. For an n × (n − k) parity check matrix H, a set of rows τ ⊂ [n] is called a stopping set if no column in H τ has one non-zero element.
Here H τ is the submatrix of H that only consists of the rows in τ .
Correspondingly, there is no parity equation that includes exactly one coded symbol among those indexed by τ . Thus, if this set of coded symbols are hidden, there is no degree-1 parity equation to recover them. Since the peeling decoder is essential for us to construct small incorrect-coding proof, the undecodable ratio α of a block is equivalent to the stopping ratio of H, which is the size of the smallest stopping set divided by n.
While CMT admits any erasure codes, SPAR uses the method introduced in [11,22] and analyzed in [26] for its proved ability to create, with high probability, parity matrices with a large stopping ratio. Given two integers c and d that satisfies nc = (n − k)d, we first generate an nc × (n − k)d permutation matrix J (a random row permutation of an nc × (n − k)d identity matrix). We then partition J into n × (n − k) slices, where each slice is a c × d sub-matrix. Then H i,j = 1 if and only if slice-(i, j) contains an odd number of 1s, for i ∈ [1 : n] and j ∈ [1 : n − k]. Such a random H has the following three critical properties: 1. It has a maximum row weight of c, and a maximum column weight of d; 2. It has a non-zero probability to have a stopping ratio of at least α * , where α * is a critical stopping ratio inherent to this method and is independent of k; 3. It is NP-hard to find the minimum stopping set and determine the stopping ratio of H. Property 1 implies that the corresponding LDPC code has a maximum parity equation size of d. Property 2 implies that we could provide the same undecodable ratio (thus same sampling requirements) for all layers. Both are desirable.
Nevertheless, Property 2 and 3 together imply that we, as the developers, are not able to determine whether the LDPC codes we generate are good (α α * ) or not, for any reasonably large k (e.g., k = 1024).
Fortunately, this problem can be easily solved through a bad-code proof. If an honest full node cannot fully decode the k data symbols after receiving (1 − α * )n coded symbols, then this code is bad, and its small undecodable set has been found and hidden by a (very strong) malicious block producer. In this case, the honest full node can prove this code bad by broadcasting the indices of the α * n coded symbols it is missing. Upon receiving and verify this bad-code proof, all the nodes in the system reject the associated block, and regenerate a code for the failed layer using an agreed random seed. This seed can be drawn from a pre-defined number sequence or the block header of a previous block, so that no consensus protocol is needed. Alternatively, distributed random number generation algorithms such as [32] and [?] can also be used for strictly unbiased randomness and consensus.
In other words, we solve the NP-hard problem of finding good codes by exploiting the computational resources of malicious party. Once it finds a small undecodable set and hides it, the system can easily detect this, reject the block, and update the code. This way, the system will settle at a good code for every layer eventually. As we will show in the next section, the probability of generating good codes is extremely high, so that SPAR can settle at good codes very quickly without having light nodes accept any fraud blocks. In addition, since bad code is a rare event, a heavy incentive/punishment scheme can be applied to deter false bad-code proof. Thus, the download and verification cost of bad-code proof is amortized to negligible throughout the course of the system.

Security
Theorem 1. In SPAR, a block producer cannot cause the soundness and agreement to fail with a probability lower than Here n i and α i are the number of symbols and undecodable ratio on the ith layer of CMT, α min min i α i , and s is the number of coded symbols each light node samples from the base layer.
Proof. Soundness: Soundness fails if a light node thinks that a block is available, and no full node is able to reconstruct the entire coded Merkle tree. We note that the reconstruction fails if any layer of the CMT cannot be recovered correctly. Let us focus on a single layer i with n i coded symbols and an undecodable ratio of α i , and assume that the malicious block producer hides α fraction of the coded symbols (and does not respond to requests for those symbols). Case-1: Consider the case of α ≥ α i . The probability of soundness failure for a node is given by the probability that a node receives all s symbols that it samples, this probability is Case-2: Consider the case of α < α i . The soundness failure occurs if a full node cannot decode the entire block or is unable to furnish a incorrect-coding proof. The full node will fail to accomplish these tasks only when it is able to receive fewer than 1 − α i fraction of symbols. Define Z i to be the total number of distinct symbols collected by the honest full node (Z i ∈ {0, 1, .., n i }). Let m be the total number of light nodes, then m · s is the total number of i.i.d. samples.

Now we have
Here (3) is by counting the number of sampling instances that provide less then (1 − α i )n i distinct symbols. H(p) = p log 1 p + (1 − p) log 1 1−p is the binary entropy function. It is apparent that we would need m large to make the above bound vanish exponentially with s.
The probability of soundness failure is smaller than the maximum probability of the two cases.
Agreement: We will argue here that soundness implies agreement for our protocol. As defined, soundness ensures that a honest full node is able to decode the block. Once a honest full node decodes the block, it will let all light nodes know that it has that block. The light nodes have either already accepted the block or have "pend"-ed the block (the light nodes could not have rejected the block since it is a valid block). The light nodes that pended the block will query with the honest full node and eventually accept the block. Thus soundness implies agreement (since now every light node agrees on the availability of the block).
Theorem 1 implies that the security of SPAR increases exponentially with the number of samples each light node takes (s), when the number of light nodes (m) is linear with the block size.

Costs and Complexity
A light node has three download costs: 1) the header, which is the CMT root of size t; 2) the random sampling, and 3) the incorrect-coding proof. In CMT, the header is the CMT root of size t. The sampling cost can be computed using the average parity-symbol-sampled Merkle proof size given in (2) to be: where b is the size of a block, and the equations hold due to that 1) s is a constant; and 2) b/k is the base layer symbol size, which is a constant. The incorrect-coding proof size can be similarly computed as where the first term is the size of d−1 coded symbols, and the second term is the size of d Merkle proofs. Finally, since the hash-aware peeling decoder decodes one coded symbol using d−1 coded symbols in one parity equation, the decoding complexity is O(1) per symbol and, thus, is O(b) in total.

Choice of parameters
Our first key parameter is the coding rate r. A smaller r means more parity symbols and thus a potentially larger undecodable ratio and less sampling. But it will also increase the height of CMT, the Merkle proof size, and the decoding complexity. For a reasonable tradeoff, we choose r = 0.25. Given r, the next two parameters we should decide are a pair (c, d) that satisfies c/d = 1 − r = 0.75 for the random LDPC code generator, where d is the maximum parity equation size. This gives us c = 0.75d and requires us to find the critical undecodable ratio of the ensemble as a function of d, which is provided in Table 2 based on the analysis in [26]. Evidently, d = 8 maximizes the critical undecodable ratio. In addition, it also admits a small incorrect coding proof that only requires only 7 coded symbols and 8 Merkle proofs. As a result, we choose (c, d) = (6, 8).

How quickly does SPAR settle at a good erasure code?
Due to random code generation, each layer of SPAR eventually settles at a good code (with an undecodable ratio of at least α * ) after a few bad codes have been deployed by the system and potentially utilized by malicious block producer to hide the data. We study the number of such attacks (note that they will never succeed) and updates before SPAR can settle: intuitively, this number can be computed as (1 − P (α < α * )) −1 − 1, where P (α < α * ) is the probability that a randomly generated code has an undecodable ratio smaller than α * = 0.124. Using an upper bound on P (α < α * ) characterized in [26] we can derive the settlement speed of SPAR as below. We note that most of the layers of CMT will immediately settle at a good code upon launching. The only exception is the layer with n = 256, which will settle after 7.7 bad codes, but without any fraudulent blocks been accepted. The proof is in Appendix A.
Theorem 2. Using the random (6, 8)-LDPC code in Section 4.4, the expected number of bad erasure codes (α < α * = 0.124) a CMT layer with n coded symbols will use before it settles at a good code (α ≥ 0.124) is approximated in Table 3.

Implementation for Bitcoin and Experiments
We developed in Rust a Coded Merkle Tree library [2] for Parity Bitcoin [33] clients (see Appendix B for more details). Our library integrates seamlessly into the current mainnet implementation, and requires minimal change on the block data structure (only need to add CMT root to the block header). Note, however, that this change is incompatible with existing Bitcoin clients. Developing a SPAR-protected Bitcoin testnet and a Bitcoin Improvements proposal (for Bitcoin Core) are ongoing research activities outside the scope of this paper.
We combine the CMT library with the performance analysis in Section 5, and numerically evaluate SPAR's light node download costs (header, sampling, and incorrect-coding proof) and full node decoding speed, for a wide range of block sizes (Table 4), and compare them with the 2D-RS based solution proposed in [9] using its C++/Python implementation [12,6]. Header ( Fig. 4(a)): A SPAR light node only downloads fixed t = 256 hashes in header, whilst 2D-RS requires 1 + 2 √ n. Thus, the header download cost of SPAR becomes much smaller than 2D-RS with growing block size. For a 64MB block, the cost is only 0.01% of the block size in SPAR, but is 0.1% in 2D-RS.
Incorrect-coding proof (Fig. 4(b)): A SPAR incorrect-coding proof only involves d − 1 = 7 coded symbols and their Merkle proofs, whilst 2D-RS requires √ k. Thus, the incorrect-coding proof download cost of SPAR becomes much smaller than 2D-RS with growing block size. For a 64MB block, the cost is only 0.051% of the block size in SPAR, but is 0.48% in 2D-RS.
Sampling cost (Fig. 4(c)): 2D-RS has a higher undecodable ratio of 25% compared to SPAR's 12.4%. Thus, for 99% confidence, s = 17 distinctive samples are enough in 2D-RS, whilst SPAR requires s = 35 if the adversary is strong enough to find, with NP-hardness, the size-0.124n stopping set. But under a realistically weak adversary that randomly selects CMT symbols to hide, SPAR only requires s = 8 because our LDPC ensemble can tolerate an average of 47% missing symbols. On the other hand, the over-sampling of each layer increases SPAR's sampling cost. Thus, although both techniques offer O(log k) sampling costs that quickly reduces with growing block size, the cost of SPAR is about 10∼16 (resp. 2.5∼4) times of 2D-RS under strong (resp. weak) adversaries. However, in practice, one can further reduce SPAR sampling cost by increasing the header size t, thus reducing the size of the Merkle proof of each symbol. Decoding speed (Fig. 4(d)): SPAR's sparse and binary encoding, at its current implementation level is already over 10 times faster than 2D-RS for all the tested block sizes.

Conclusion and Discussions
By iteratively applying a special ensemble of LDPC codes to every layer of a Merkle tree and batching the hashes of each coded layer into the data symbols of the next layer, we invented a novel hash accumulator called coded Merkle tree (CMT). Built upon CMT, we proposed a novel data availability verification system called SPAR, which allows the availability and integrity of the entire tree to be checked at constant costs.
SPAR can play a key role in scaling blockchain systems that incorporate light nodes because it empowers these nodes with real-time verification of data availability and integrity at small and constant costs. SPAR can also be used to scale the communication of sharded blockchain systems (e.g., [23,19,20]), where full nodes of one shard operate as light nodes of other shards, as SPAR allows them to efficiently check the availability and integrity of blocks in other shards.
Integrating SPAR into existing blockchain systems requires minimum changes and no extra bandwidth consumption. An honest block producer only needs to broadcast the original data block as usual and attach the CMT root in the block header. This is sufficient for other full nodes to reproduce the CMT and offer sampling services for light nodes. Our library for CMT in Rust for Parity Bitcoin clients maintains the same API as the standard Merkle tree module. Noting that classic Merkle trees are indeed special CMTs with coding rate r = 1 and batching factor q = 2, our library readily replaces the standard module and is backward compatible. the probability that the stopping distance of the code is smaller than α * n is upper-bounded by P (αn < α * n) ≤ min inf 0<δ<α * n(α * − δ)e n max θ∈[δ,α * ] γ(θ) where for c = d(1 − r), h(θ) = −θ log θ − (1 − θ) log(1 − θ), and x 0 as the only positive solution to For α * = 0.124, the above upper bound for small n (e.g., n = 256) becomes degenerated (i.e., reduces to the trivial bound of 1). In order to obtain a good approximation for all the values of n, we approximate the upper bounds on P (α < 0.124) using a slightly smaller undecodable ratio of 0.116. Then, we evaluate the upper bounds on P (αn < 0.116n) in (8), for all the considered values of n in Theorem 2 to obtain the probabilities in the second row of Table 3.
We note that since 0.116 is very close to 0.124, SPAR's inherent oversampling of intermediate layers will provide sufficient protection for data availability on these layers, so that the light node sampling cost will not increase.

B Coded Merkle Tree Library
We developed in Rust a Coded Merkle Tree library [2] for Parity Bitcoin [33] clients. We modify the data structure of the block header to add a new field coded merkle roots hashes, which are the hashes of the coded symbols on the last level of the coded Merkle tree constructed from this block.
To use the Coded Merkle Tree library on a block, we require the following input parameters from the users: -BASE SYMBOL SIZE: size of a symbol on the base level, and hence the number of systematic symbols on the base level. -AGGREGATE: number of hashes to aggregate into a symbol on the next level.
-HEADER SIZE: number of hashes stored in the block header. This also decides the total number of levels in the coded Merkle tree. -Codes for all levels of coded Merkle tree, in forms of sparse representations of their parity-check matrices. Given the above parameters, Coded Merkle Tree implements the following key functionalities: coded merkle roots: construction of the coded Merkle tree from the block content.
merkle proof: generating the Merkle proof for any symbol on any level of coded Merkle tree. By design, this returns a set of symbols on the higher level. sampling to decode: sampling symbols on the base level, together with their Merkle proofs. run tree decoder: decode the entire coded Merkle tree level by level from the roots. Each level is decoded by running a hash-aware peeling decoder, using the decoded symbols on the previous level as the hash commitments. generate incorrect coding proof: 1) when a coding error is detected, this function returns d − 1 symbols in a parity equation, and Merkle proofs for all d symbols in that equation; 2) when the peeling process gets stuck before all symbols are decoded, this function returns the indices of the missing symbols as a stopping set.


extfieldsofstudy:
Mathematics
Computer Science

text:
Total Roto-Translational Variation

We consider curvature depending variational models for image regularization, such as Euler's elastica. These models are known to provide strong priors for the continuity of edges and hence have important applications in shape-and image processing. We consider a lifted convex representation of these models in the roto-translation space: In this space, curvature depending variational energies are represented by means of a convex functional defined on divergence free vector fields. The line energies are then easily extended to any scalar function. It yields a natural generalization of the total variation to the roto-translation space. As our main result, we show that the proposed convex representation is tight for characteristic functions of smooth shapes. We also discuss cases where this representation fails. For numerical solution, we propose a staggered grid discretization based on an averaged Raviart-Thomas finite elements approximation. This discretization is consistent, up to minor details, with the underlying continuous model. The resulting non-smooth convex optimization problem is solved using a first-order primal-dual algorithm. We illustrate the results of our numerical algorithm on various problems from shape-and image processing.

Introduction
It was observed at least since [58] that line energies such as Euler's "Elastica" (here for a smooth curve γ): γ [1+] κ 2 ds could be natural regularizers for the completion of missing contours in images. This idea was based in particular on observations of G. Kanizsa [46,47] about the way our perception can "invent" apparent contours. The model in [58] (see also [40,39] for interesting attempts to solve it with phase-field methods) was variational in nature and since then, many attempts have been proposed to study and address the minimization of such line energies, both theoretically and numerically.
From the theoretical point of view, the lower-semicontinuity of these energies already is a challenge, which has been studied in many papers (in connexion to the applications to computer Riemannian and sub-Finslerian geodesic tracking in the roto-translation group [11,34]. As said, we wish to extend these variational methods to level sets representations. The family of problems we are interested in is introduced in a paper by Masnou and Morel [52] (see also [53]), which addresses the problem of image inpainting. Their initial idea is to minimize, in the inpainting domain D, an energy of the form D

+ div
Du |Du| p |Du| for some p ≥ 1, where u is a bounded variation function and appropriate boundary conditions are given. A theoretical study of this energy, which in general is not lower-semicontinuous (lsc), is found in [4], in particular it is shown that if p > 1 (in dimension 2), it is lsc on C 2 functions. Interestingly, also this study relies on a space/direction representation (and more precisely on varifolds). An interesting co-area formula for the relaxed envelope is also shown in [51].
There have been many attempts to numerically solve Masnou and Morel's model, which is very difficult to tackle, being highly non convex. Most of the authors introduce auxiliary variables [9,26,10], for instance representing the orientation, which is already close to the idea of the Gauss map or Roto-translational representation. Recent techniques based on Augmented Lagrangian methods (for coupling the auxiliary variables) [74,80,78,79,43,42,8] have shown to be quite efficient, despite the lack of convexity and hence convergence guaranties.
What we propose here is to rather introduce a functional which may be sees as a convex relaxation of Masnou and Morel's. The method we propose is based on two ingredients: A lifting in the Roto-translation space of curves which allows to write the Elastica energy (or any convex function of the curvature) as a convex function, as in classical Gauss map based approaches, and (formally) a decomposition of functions as sum of characteristic functions of sets with finite energy. This allows to define a convex functional which is defined on grey-level functions and penalises the curvature of the level lines. It is however easy to check that it is in general strictly below Masnou and Morel's functional, in particular a function of finite energy needs not have its level sets regular (in practice, they could be intersections of regular sets, in addition to having possible cusps, for the same reasons as in [12]), see for example Figures 5 and 9(e). The only theoretical result which we can show is that C 2 curves are tightly represented by our convex relaxation.
There is a close relationship between our approach and the functionals in [22] ("TVX") and [23], based on similar representations (but [22] omits to preserve the boundary of the lifted current). In fact, the functional we build is, as we show further on, a new expression of the previous convex relaxation of the Elastica energy (and variants) proposed in [23]. This stems from the identity of the corresponding dual problems. Although it boils down to the same energy, it is introduced in a much simpler way, as the primal expression in [23] requires to work in a space where the point, the tangent and the curvature are lifted as independent variables. As a consequence, we also can provide a simpler discretization, and the tightness of the relaxation for C 2 sets was unnoticed in [23]. Eventually, we should mention that part of what we propose here could be generalized to arbitrary dimension or co-dimension, following the techniques in [6]. However, computationally, our construction for 1D curves in the plane already needs to work in a 3D space, and more complex models seem at this moment intractable. Our future research will rather focus on improving the discretization and the optimization of the bidimensional case.
The paper is organised as follows: In Section 2 we describe the lifting and introduce our functional. We state our main result, which shows that characteristic functions of C 2 sets are well represented by our convex relaxation (Theorem 1). In the next Section 3 we show an approximation result with smooth functions, and compute a dual representation of the functional (which shows it is identical to [23]). Then, in Section 4, we describe how the functional is discretized and show a few experimental results. The last Section 5 is devoted to the proof if Theorem 1. Some technical tools are found in the Appendix, in particular, Appendix A shows, up to (hopefully minor) transformations, the consistency of our implementation in Section 4.

The functional
We consider Ω a bounded open domain in the plane, and E ⊂ Ω a set with C 2 (and to simplify, connected) boundary. Assume that this boundary ∂E is given by a parameterized curve (x 1 (t), x 2 (t)) with parameter t ∈ [0, 1]. The (extrinsic) curvature κ E of ∂E is classically defined as the ratio between the variation of the tangential angle θ and the variation of its arc length s, that is The main idea of the lifting is now to consider a higher dimensional representation of the parametric curve in the 3D roto-translation (RT) space, which is obtained by adding the tangential angle θ as an additional dimension to Ω. In this space, we now consider a parameterized 3D curve (x 1 (t), x 2 (t), θ(t)) which lifts the boundary ∂E to the roto-translational space. Figure 1 shows an example where we lift the boundary of a disk (which is a circle) to the RT space. Observe that the lifted boundary is represented by a 3D helix. Now, we define for all t ∈ [0, 1] the tangential vector p(t) = (p x (t), p θ (t)) with The curvature is therefore given by In this work we consider f : R → [0, +∞] a convex, lsc function and want to define a convex lsc extension (to grey-level valued functions) of energies of the type where E ⊂ Ω is a set with C 2 boundary, and κ E is the curvature of the set. Using our tangential vector p(t), it is easy to see that the energy can be (formally) written as 1]) is the lifted curve and τ (x, θ) its normalized tangential vector, given by τ (x(t), θ(t)) = p(t)/|p(t)| for all t ∈ [0, 1]. A precise definition of the energy will be given below.
Let us briefly discuss three instances of energies that will typically appear in applications. In all cases α > 0 will be a tuning parameter that can be used to balance the influence of the curvature term with respect to the length. Figure 1: The gray line is the lifting of the boundary ∂E of the disk E to the roto-translational space Ω × S 1 .
1. f 1 (t) = 1 + α|t|. This energy penalizes the arclength plus the absolute curvature, hence we might expect that this type of energy will allow also for corners. This type of energy has been studied with a different approach in [22]. The energy is given by: 2. f 2 (t) = 1 + α 2 |t| 2 . This energy penalizes the arclength of the lifted curve in the RT space which, in some sense, corresponds to the "Total Roto-translational Variation". It yields the energy: which is nothing but the length of the lifted curve Γ E in a Riemanian metric. This functional was considered for perceptional completion problems probably first in [65]. It coincides with sub-Riemannian [65,33,20,11] and more precisely the sub-Finslerian models in [34] (with the constraint of positive direction on the velocity).
Observe that the second term in the energy is a "quadratic over linear" function which is still convex, provided τ x is constrained in a half space or a half line (which of course will depend on the θ variable).
In all these examples, we find that the energy which we are interested in ends up represented possibly as a convex function of the measure σ = τ H 1 Γ E . Moreover, this σ is not an arbitrary measure. In particular, it satisfies two important constraints: first, by construction, it is a circulation and has zero divergence in Ω × S 1 (it can have source terms on ∂Ω × S 1 ): indeed, given any smooth function ψ, ∇ψ · σ = Γ E ∂ τ ψdH 1 vanishes if Γ E is a closed curve, or if it ends on the boundary and ψ has compact support. Second, its marginals in Ω × S 1 , which we can formally denote S 1 σ (and which are also divergence free), coincide with a 90 • -rotation of the measure Dχ E . We will show now how to generalize this construction to arbitrary sets or functions (with bounded variation).
In the whole paper, we assume that there exists γ > 0 such that We also introduce the recession function which is a convex, one-homogeneous function possibly infinite on (−∞, 0) or/and (0, +∞). It is easy to check that it is the support function of dom f * , the domain of the convex conjugate f * of f : Indeed, for y ∈ dom f * , s > 0, as s → ∞; on the other hand, as f * ≥ −γ thanks to (1), for s > 0 one has We let then, for p = (p x , p θ ) ∈ R 3 with p x = 0, Here, θ denotes the unit planar vector (cos θ, sin θ) T (by a slight abuse of notation, we will also denote in this way the vector (cos θ, sin θ, 0) T ∈ R 3 ). It is then classical [64, § 13] and easily follows from (2) that one has h(θ, p) = sup ξ · p : that is, the one-homogeneous function h is the support function of the convex set in the right-hand side of (4). Indeed, first of all, the sup in (4) is ∞ if p x is not λθ, λ ≥ 0. While if it is of this form, then (taking the suprema over the ξ's which satisfy ξ Observe that (1) yields that for all p ∈ R 3 , We now introduce the functional Here x ⊥ = (x 2 , −x 1 ) is a 90 • rotation in the plane. The last condition is understood as follows: In particular in general the fields σ appearing in (6) are free-divergence bounded Radon measures in Ω × S 1 with values in R 3 (we denote M 1 (Ω × S 1 ; R 3 ) the space of such measures) and the proper way to write the integral is rather where σ/|σ| is the Radon-Besicovitch derivative of σ w.r. its total variation |σ|. Notice that for any u and any admissible σ, thanks to (5), because of (7), taking the supremum with respect to the test function ϕ with |ϕ(x)| ≤ 1 everywhere. This shows that F (u) bounds the BV seminorm. We denote by K(θ) the closed convex set of R 3 whose support function is h(θ, ·), as already seen it is given by Remark that (1) implies that f * (s) ≤ − γ 2 − s 2 for |s| ≤ γ, implying in particular that 0 is in the interior of K(θ). We then let (following for instance [21,63]) Then, we claim that for any measure σ ∈ M 1 (Ω × S 1 ; R 3 ) showing in particular that this integral is a lower-semicontinuous function of σ (for the weak- * convergence). To prove this claim (which is standard), first observe that if ϕ ∈ K, then for all ψ − ∈ C 0 c (Ω × S 1 ; R + ) and ψ ∈ C 0 c (Ω × S 1 ), also ϕ − ψ − θ + ψθ ⊥ ∈ K. One deduces easily that the supremum in (8) is infinite if one of the measures (θ · σ x ) − or θ ⊥ · σ x does not vanish. If both vanish, it means that σ = (λθ, σ θ ) for some nonnegative measure λ. Hence the supremum becomes and (8) is then deduced in a standard way (with a Besicovitch covering of Ω × S 1 with respect to the measure |σ| by balls where λ, σ θ are "almost constant", and choosing then for ϕ the "right function" in each ball of the covering-the fact that zero is in the interior of K, which allows to multiply functions in K by a cut-off, is important here). We can deduce that also F is lower semicontinuous: consider a sequence (u n ) with F (u n ) ≤ c < ∞. Then u n is bounded in BV and converges (up to a constant and a subsequence) to some u in L 1 (Ω).
If σ n reaches the value of F (u n ) up to 1/n, one has that Ω×S 1 |σ n | ≤ c < ∞ so that σ n * σ (up to a subsequence), as measures, and obviously div σ = 0. Clearly, also (7) passes to the limit. Hence, by lower semicontinuity, This shows that F defines a convex, lower semicontinuous functional on BV (Ω). Our main theoretical result is the following theorem, which shows that if the argument u is the characteristic function of a smooth enough set, then F (u) coincides, as expected, with a curvature-dependent energy of the boundary of the set {u = 1}. The proof of this result is postponed to Section 5.
Remark 2.1. One could hope that F coincides with the lower semicontinuous envelope of its restriction to C 2 sets, with respect to the L 1 convergence. However, simple examples show that it is not the case. Many examples where it fails are found in [23]. In general, we expect the relaxation F to be strictly below the relaxations found in the literature since [12], based on the (L 1 ) approximation of sets with smooths sets. In particular, in the case of [12,Fig. 4.2], our relaxation will certainly be below twice the energy reached by the pattern in [12, Fig. 1 2. If f has growth one and is even, one expects that the same result holds for piecewise C 2 sets. Remark 2.3. We believe that the proof below could be extended to show that if u ∈ BV (Ω) is a function with C 2 level sets, and such that there exists a continuous function κ(x) which coincides with the curvature of ∂{u > s} for all level s ∈ R, then one should have Hence in that case our functional would coincide with Masnou and Morel's [53]. A typical example of a function u for which such a function κ(x) exists is (in a convex domain Ω, see [54]) a solution of the so-called "Rudin-Osher-Fatemi" functional for some λ > 0, in case f is continuous. Then, also u is continuous and κ(x) is given by λ(u−f ). In that case, in addition, κ(x) is also the curvature of the level sets of u = h(u) for any nondecreasing function h such that h(u) is still in BV (Ω).
3 Some properties of the functional 3.1 Approximation by smooth functions Assume Ω is a bounded convex set. Then for any u ∈ L 1 (Ω) with F (u) < ∞, there exists (u n ) n a sequence of functions with u n ∈ C ∞ (Ω), which converge to u in L 1 (Ω), and such that lim Remark 3.2. We believe that, upon replacing C ∞ (Ω) with C ∞ (Ω), the result should be true in any domain.

Dual representation
We can show the following representation for F , which in particular establishes that F coincides exactly with the (more complicated) relaxation R * * el defined in [23].
Assume ∂Ω is connected. Then the functional (6) can also be represented in the following form If ∂Ω is not connected, then same formula holds, however the functions ϕ should be in C 1 (Ω × S 1 ), with ∇ϕ ∈ C 0 c (Ω × S 1 ).
Corollary 3.4. The functional F (u) coincides with the relaxation of the Elastica energy R * * el (u) defined in [23].
Indeed, formula (12) is the same as (35) in [23,Theorem 13]. The primal formulation (6) we introduce here is quite simpler as the original formulation in [23], though, as the "κ" variable is implicit in our formulation (it comes naturally as the derivative of the orientation) and does not need to be lifted.
Proof of Prop. 3.3. We use the standard perturbation approach to duality [35], in the duality (C 0 0 (Ω × S 1 ; R 3 ), M 1 (Ω × S 1 ; R 3 )), where C 0 0 denotes the L ∞ -closure of the functions with compact support (hence, the functions which vanish on the boundary) and M 1 the totally bounded vector-valued Radon measures. Defining, for p ∈ M 1 (Ω × S 1 ; R 3 ), the function so that F (u) = G(0), we observe that G is lower semicontinuous. Indeed if p n * p and σ n is a minimizer in the definition (13) of G(p n ), as before thanks to (5), σ n + p n is bounded and up to a subsequence, converges to a measure σ + p where σ is admissible in (13). It follows that In particular, we deduce that G = G * * so that It remains to compute G * (η), to show that (12) follows: and we find that Recalling that the sup is over the σ admissible, we find that Now, given any σ with div σ = 0, we can define u such that hence the first supremum in (14) is the same as the sup over all σ with vanishing divergence. Hence, it is zero or +∞, depending on whether η − ψ is a gradient or not. We find eventually G * (η) is finite only when there exists ϕ ∈ C 1 (Ω × S 1 ) with ∇ϕ = 0 on ∂Ω × S 1 , such that η = ψ + ∇ϕ. The proposition follows.

Numerical Experiments
For numerical solution of the proposed model, we need to discretize both the 2D image domain as well as the 3D domain of the roto-translation space. Due to the high anisotropy of the energy, one has to be extremely careful in the choice of the discretization scheme in order to preserve a maximal degree of rotational invariance while keeping the numerical diffusion as low as possible. It turns out that a 2D-3D staggered grid version of an averaged first-order Raviart-Thomas divergence conforming discretization [62] yields the best results. More elaborate discretizations using for example adaptive grids or higher-order approximations will be subject for future study. Figure 2 shows a visualization of our combined 2D-3D staggered grid approach. The first grid is given by a 2D grid of pixels discretizing the image domain Ω. The second grid is given by a 3D The bottom grid represents the 2D image, and the image intensities u are stored at the vertices of the squares S i,j . The gray arrows indicate differences between the pixels. The volume V i,j,k represents one element of the corresponding discretization of the RT space. The black arrows through the faces of the volume represent the 3D vector field σ.

Staggered averaged Raviart-Thomas discretization
grid of volumes which discretize the roto-translation space Ω × S 1 . In the visualization we show a few image pixels of the 2D grid (gray nodes) and one cube of the corresponding 3D grid. First, we start by describing the discretization of the image domain. We assume here that Ω is a square or rectangle, or in general a convex domain, which we will discretize at a scale δ x > 0. We will restrict ourselves to the case of a rectangle to simplify the notation. Then, we let Ω δx denote the discretized image domain which is given by denotes the square centered around the index i = (i, j). The intensity values of the discrete image are stored at the vertices half-grid points (such as the four vertices of S i,j ), and for notational simplicity we introduce the following index set: Hence, the discrete image u ∈ R I consists of N 1 × N 2 discrete pixels with values (u i ) i∈I . We will later also need the following index sets: which index the middle points of the edges between two adjacent pixels along the first and second spatial dimensions. Next, we let δ θ > 0 be the discretization step for the angular variable. We must assume that We then denote by Γ δ with δ = (δ x , δ θ ) the discretized 3D roto-translation space. It is given by denotes the volume centered around the 3D index j = (i, j, k). We associate an index set J of center locations of the volumes V i,j,k : We shall also introduce the following three index sets: which we will use to index the facets of the volumes. Observe that the 3D points which are orthogonal to the first spatial dimension. Similarly, the 3D points (iδ x , (j − 1 2 )δ x , kδ θ ) for all (i, j − 1 2 , k) ∈ J 2 correspond to the middle points of the facets which are orthogonal to the second spatial dimension. Finally, the 3D points (iδ x , jδ x , (k − 1 2 )δ θ ) for all (i, j, k − 1 2 ) ∈ J θ correspond to the middle points of the facets which are orthogonal to the angular dimension. Now, we consider a discrete vector field σ = (σ 1 , respectively. The discrete vector field σ can also be seen as a Raviart-Thomas vector field [62], defined everywhere in Γ δx,δ θ . It is obtained by an affine extension: where ∆(t) = max(0, 1 − |t|) is the usual linear interpolation kernel.
It is well know that such a Raviart-Thomas field has a divergence (defined everywhere in the distributional sense) which is obtained in each volume by summing the finite differences of the average fluxes through each pair of opposite facets. In our model, the discrete field σ is constrained to be divergence free, hence for each volume V j , we must impose that: The roto-translation space is periodic in the θ direction and hence, whenever k = N θ , we identify the index N θ + 1 2 with the index 1 2 in the computation of the last term of the discrete divergence. To simplify our presentation, we introduce a discrete divergence operator D : R J 1 ∪J 1 ∪J θ → R J , so that the divergence free constraint (22) can be compactly written as Next, to implement a discrete version of the consistency condition we observe that averaging (in the direction of θ) the first two components of the field σ yields a rotated (by π/2) gradient which should locally coincide with the finite differences of the discrete image u. This yields the following discrete consistency condition: We introduce a projection operator P : R J 1 ∪J 2 → R I 1 ∪I 2 and a discrete (rotated) gradient operator G : R I → R I 1 ∪I 2 such that we can write the above discrete compatibility condition as In order to approximate the continuous energy (6) with our discrete Raviart-Thomas field σ RT we can use different types of quadrature rules. After various attempts, we found out that simply summing the values at the center points of the volumes V i,j,k provides the highest flexibility for the discrete field σ to concentrate on thin lines, and, in turn, yields the most faithful numerical results.
Letting thus σ i,j,k = σ RT (iδ x , jδ x , kδ θ ), we have that this discrete field is obtained by the following formula: ). (24) We again introduce an operator A : R J 1 ∪J 2 ∪J θ → (R 3 ) J such that the above averaging operation can be writtenσ = Aσ.
Now, using the volume centered quadrature rule, the discrete energy is The consistency of this energy, as δ := (δ x , δ θ ) → 0, with the continuous energy F defined in (6) is studied in Appendix A.
The regularization function h(θ, p), where θ ∈ S 1 and p = (p x , p θ ) ∈ R 3 is defined as in (3), and the function f (t), which appears in its definition will be one of the following classical examples of convex functions: In all examples, the regularizing function is a combination of length and curvature regularization and the parameter α > 0 can be used to adjust the influence of the curvature regularization.
The first function, f 1 is the sum of length and absolute curvature, hence we call the corresponding regularizer "total absolute curvature" (TAC). One of its main features is that it allows for sharp corners in the level sets of the image. Interestingly, since integrating the absolute curvature along the boundary of a shape is constantly 2π for all convex shapes (and in general, is scale independent) then the main effect of the curvature term in this energy is to penalize non-convex shapes, regardless of their size.
The function f 2 combines length and curvature through an Euclidean metric and hence, it corresponds to the total variation of the lifted curve in the RT space, hence we consequently denote this regularizer "total roto-translational variation" (TRV). For relatively small curvature, it favors smooth shapes, but it also allows sharp discontinuities.
Function f 3 penalizes squared curvature plus length and hence is equivalent to the Elastica energy. We call this regularizer "total squared curvature" (TSC). It is well-known that while length regularization favors smaller shapes, quadratic curvature favors larger shapes. Hence, the interplay between length and curvature regularization removes the shrinkage bias and leads to smooth shapes.
Based on the "basis" functions f l , l = 1, 2, 3, we obtain the following convex regularization functions in the RT space.
Observe in particular, that the above definitions already properly take into account the correct values of the function h l , l = 1, 2, 3 for p x = 0. In the following sections we will apply the proposed curvature based regularization functions to a variety of image-and shape processing problems. For this, we consider generic optimization problems of the form min where G(u) is a convex, lsc. function defining an image-based convex data fidelity term, possibly dependent on the pixel location.

Primal-dual optimization
In this section, we show how to compute a minimizer of the non-smooth convex optimization problem (29). For notational simplicity, we first divide the whole objective function (29) by δ 2 x δ θ and assume that the remaining factors are absorbed by the data fitting term.
In order to solve the constrained optimization problem (29), we consider its Lagrangian (saddlepoint) formulation: where φ ∈ R J , ξ = (ξ 1 , ξ 2 , ξ θ ) ∈ (R 3 ) J and ψ = (ψ 1 , ψ 2 ) ∈ R I 1 ∪I 2 with ψ 1 ∈ R I 1 and ψ 2 ∈ R I 2 are the dual variables (Lagrange multipliers or discrete test functions). The function h * denotes the convex conjugate of the function h. Recall that h is the support function of the convex set where f * denotes the convex conjugate of f . Hence, the convex conjugate h * is simply the indicator function of the set H: The problem (30) is a saddle-point problem which is separable in the non-linear terms and hence falls into the class of problems that can be solved by the first-order-primal-dual algorithm with diagonal preconditioning and overrelaxation [24,60,25]. In order to make the algorithm implementable, we need efficient algorithms to compute the projection operators proj H(θ) and proximity operators prox τ G . They are defined as the unique minimizers of the minimization problems: proj H(θ) (η) = arg min and, for some τ > 0. In what follows, we will detail the projection operator for different instances of the convex sets H (respectively regularization functions h), the proximity maps for G will be detailed as soon as they are needed in the numerical results. The general idea for performing the projection ξ = proj H(θ) (η) of a point η = (η x , η) onto the set H(θ) is outlined in Figure 3. The Figure on the left hand side represents the set H in the ξ x plane, and the Figure on the right hand side shows the "profile" of the set which is obtained by cutting along the ξ x θ × ξ θ plane. We let η x θ = η x · θ and define (ξ x θ , ξ θ ) = proj P (η x θ , η θ ) as the projection of the point (η x θ , η θ ) onto the "profile" Then, using simple geometric reasoning, we see that the variable ξ x can be recovered from ξ θ by It remains to detail the projections onto different profiles P . Figure 3: Projection of the point (η x , η θ ) onto the convex sets H. First, the point (η x θ , η θ ) is projected onto the profile P which yields (ξ x θ , ξ θ ), then ξ x is computed from ξ x θ via the geometric relationship ξ • TAC: f 1 (t) = 1 + α|t|: The convex conjugate of the function f 1 is given by and in turn the profile P 1 is given by It is straightforward that the projection of a point (η x θ , η θ ) onto the profile P 1 can be performed via simple truncation operations .
• TRV: f 2 (t) = √ 1 + α 2 t 2 : A simple computation shows that the convex conjugate of f 2 (t) is given Inserting the expression of the convex conjugate into the profile (31) and squaring both sides, we obtain In what follows, we assume that max(0, η x θ ) 2 + (η θ /α) 2 > 1 since otherwise we do not need to project the point. We first treat the case η x θ ≤ 0. It is easy to see that in this case the solution of the projected point is given by min(α, η θ )) .
In case η x θ > 0, computing the projection of a point (η x θ , η θ ) onto the boundary of P 2 amounts to solve the following equality constrained optimization problem The Karush-Kuhn-Tucker (KKT) optimality conditions for the above problem are given by where λ > 0 is the Lagrange multiplier which is positive since the point (η x θ , η θ ) was assumed to be outside P 2 . Combining the first three equations shows that the Lagrange multiplier is computed from the roots of the fourth-order polynomial Let us observe that in case α = 1, the solution for λ is particularly simple, indeed such that from the first two equations we obtain for the projected point In the general case α = 1 we compute λ by applying Newton's algorithm to find the (correct) root from the fourth order polynomial. In our experiments it turns out that a large enough initial value, e.g. λ = 10 3 provides a good initialization for Newton's algorithm. Usually, we need less than 5-10 iterations of Newton's algorithm to converge to a solution with feasibility error less than 10 −9 . From the first two equations of the KKT conditions the expression of the projected point is given by: .
• TSC: f 3 (t) = 1 + α 2 t 2 : The convex conjugate of f 3 is computed as and hence the profile P 3 is given by Following the same approach as before a point (η x θ , η θ ) with η x θ + (η θ /(2α)) 2 > 1 is projected onto P 3 by solving the KKT optimality conditions Combining the above equations, the Lagrange multiplier is obtained by finding the correct root from the third-order polynomial We again apply Newton's method and observe that for a large enough initial value, e.g. λ = 10 3 , Newton's algorithm rapidly converges to the correct root of the polynomial. The projected point is finally given by

Computing a disk
In the first example, we consider the most basic numerical experiment, which is using TSC energy to compute a disk of a given radius. The aim of this experiment is to investigate the quality of our proposed discretization scheme, in particular when using a different number of discrete orientations. Consider a 2D disk B(0, r) = {x ∈ R 2 : |x| ≤ r} of radius r > 0 centered around the origin. The TSC energy of the boundary of the disk B(0, r) is given by Here our disk will be represented by its characteristic function χ B(0,r) (x). We can force the minimizer of the TSC energy to yield a disk using suitable boundary conditions. For example, we can force at least one point inside the disk to be one and at least one point outside the disk to be zero. The following simple computation shows that the minimizer of the TSC energy will be a disk of radius r = α, indeed: which is minimized if κ = 1/r = 1/α everywhere. We set up a our optimization problem using a 2D grid of N 1 × N 2 = 40 × 40 pixels. We will use a different number of discrete orientations N θ ∈ {4, 8, 16, 32, 64} in order to investigate the quality of the approximation depending on N θ . In this and the subsequent experiments, the discretization width of the spatial grid is set to δ x = 1 and the discretization width of the angular dimension is set to δ θ = 2π/N θ . The aim of this experiment is to compute a disk of radius r = 10 pixels. We therefore define an inpainting domain D ⊆ I forming a band of 10 pixels width around a disk of radius r = 10. In this domain we minimize the TSC energy and we use the remaining part of the image as boundary condition. The corresponding data term in our optimization problem (29) is given by where ι C denotes the indicator function of the convex set C. The proximal map for G(u) is given by  In order to quantify the approximation quality of our discretization for a varying number of discrete orientations, we report the total variation (TV), the absolute curvature (AC) and the squared curvature (SC). All three quantities are computed from the averaged fieldσ = Aσ using where the function h is one of the three instances: . Table 1 details the values of the discrete energies we obtained for different numbers of discrete orientations N θ . From the results, one can see that a higher number of discrete energies generally leads to a better approximation of the true energy. This is particularly true for the value of the squared curvature which seems to be well approximated only when using a quite high number of discrete orientations. The absolute curvature, however seems to be well approximated even when using only a small number of discrete orientations. The reason for this is that the absolute curvature of a smooth curve is easy to approximate by means of a piecewise linear curve. Figure 4(f) visualizes the averaged Raviart-Thomas vector fieldσ of the disk example in case of N θ = 64 discrete orientations. Observe that the vector field nicely corresponds to the expected shape of a helix, shown in Figure 1. However, due to diffusive effects of our numerical scheme, the vector field does not perfectly concentrate on a one-dimensional structure.

Non-smooth level sets
Here we demonstrate a typical effect of our convexification which finds a low-energy solution for an image with non-smooth level sets which should have infinite energy in more standard relaxations [53] of the Elastica energy. Figure 5(a) shows an image of size N 1 × N 2 = 100 × 50 pixels of a black and a white disk in front of a gray background. Similarly to the inpainting problem of the previous example, we fix only the four small parts of the original image (see Figure 5(b)) and minimize the TSC energy in the inpainting domain D, indicated by the gray area. In this experiment we used N θ = 64 discrete orientations and we set α = 17 to match the radius of the disks. Figure 5(c) shows the computed minimizer of the TSC energy. Observe that the solution u does not yield the expected two disks of the original image but rather drop-like shapes that form a sharp cusp in the middle of the image. Figure 5(b) shows a stream-line representation of the  minimizing vector fieldσ in the roto-translation space. Inspecting the field, one can immediately see the reason for this behavior. The vector fieldσ forms a twisted ∞-shape curve which "skips" the strong curvature of the cusp.

Shape completion
In this section we provide some qualitative results on a number of different shape completion problems. Similar to the previous two examples, we define an inpainting domain D which is indicated by the gray area and keep the remaining image as boundary condition. Figure 6 shows various input shapes with their inpainting domains and the solutions of minimizing different curvature energies using different settings of the parameter α. One can see that while minimizing the TAC energy usually leads to straight connections and sharp corners, TSC leads to a smooth continuation of the boundaries. The TRV energy leads to results which are somewhere in between the results of TAC and TSC. In the last row of Figure 6, we additionally demonstrate the behavior for completing straight lines at different rotations by minimizing the TAC energy. It turns out that  : Shape completion for various shapes using total absolute curvature (TAC), total rototranslational variation (TRV) and total squared curvature (TSC). The inpainting domain is indicated by the gray pixels. In the last row we provide results for completing straight lines at different angles.
our discretization scheme performs quite well for rotations of 0, π/4, π/2, . . . but leads to more diffusive results for rotations of π/8, 3π/8, . . .. The development of a more isotropic scheme will be subject of future research. In Figure 7, we show the application to shape completion from dipoles using the original data of "Weickert's cat" [76]. The input image u 0 is of size N 1 × N 2 = 400 × 380 pixels and the inpainting domain D is again indicated by the gray pixels. In this example we used N θ = 64 discrete orientations and the curvature parameter was set to α = 50. In order to avoid relaxation artifacts (self-intersections) inside the dipoles, we additionally use zero-boundary conditions for the field σ 1,2 within the constant areas of the dipoles. The results show that our numerical scheme can successfully reconstruct the curvilinear shape of the cat from only very little information given by the dipoles. However, we can also observe diffusion artifacts of our numerical scheme, especially if the reconstructed boundaries are relatively long. Sharper results for this kind of problems are usually obtained using sophisticated anisotropic diffusion schemes such as the edge enhancing diffusion (EED) [75,67]. It would be interesting to understand whether these schemes also minimize an underlying variational energy.

Shape regularization
In our next experiment, we apply our curvature based energies for shape regularization. Given an input image (u 0 i ) i∈I , which is the characteristic function of a given shape, our aim is to compute a simplified (or regularized) shape which is represented by means of a (relaxed) binary image (u i ) i∈I . We make use of a simple linear fidelity term which is frequently used in image segmentation: where w ∈ R I is a force field. For the application to shape regularization we use w = λ( 1 2 − u 0 ), where λ > 0 defines the strength of the force field. The proximal map for this data term is easily computed: In Figure 8, we apply TAC, TRV and TSC regularization to regularize the shape of Picasso's "Bull fight" image. In all three cases, we set α = 10 and we use different settings of the parameter λ to obtain gradually simplified shapes. From the results one can see that TAC regularization yields shapes with relatively straight boundaries and sharp corners. TRV regularization yields smooth shapes but also allows for sharp corners. TSC regularization yields smooth shapes. Observe that whenever it seems energetically preferable, the solution of TSC produces sharp cusps with "hidden" edges to bypass locations of strong curvature, as predicted by the theory [12]. This effect is usually less visible when minimizing the TAC or TRV energies.

Image inpainting
In this section, we apply our proposed curvature energies to the classical problem of image inpainting. Similar to the previous examples, we use the data term (32) with the only difference that the input image u 0 is now a gray level image.
In both examples we used N θ = 32 discrete orientations and the curvature parameter was set to α = 15. Figure 9 shows the results for two different inpainting problems. In the "Einstein" image, we randomly remove 80% lines and in the "Picasso" image we randomly remove 90% of the pixels. In case of the "Einstein" image, we use TSC regularization with α = 10 and for the "Picasso" we set α = 5. In case of the "Einstein" image, the gaps a much larger and hence we use a larger value of α. For comparison, we also provide results of standard TV regularization, which is equivalent to using α = 0 in one of the three curvature energies. From the results one can clearly see that curvature regularization leads to significantly better inpainting results. On the downside, we can also observe some artifacts which are caused by our convex representation in the roto-translation space.

Image denoising
In our last experiments we investigate our curvature energies for the classical problem of image denoising. We investigate two different types of noise: Zero-mean Gaussian noise and impulse noise such as "salt & pepper" noise. For Gaussian noise, it is well-known data a quadratic 2 data term where λ > 0 is a good data fidelity parameter. The proximal map is given by In case of impulse noise, a 1 data term is more suitable: since it is more robust with respect to outliers. The proximal map for the 1 data term is given by the classical soft-shrinkage formula In Figure 10 we show the results of our curvature-based regularization energies for image denoising. In case of the "Louvre" image we generated the noisy image by adding 25% "salt & pepper" noise. The noisy "Leberblümchen" image was generated by adding zero-mean Gaussian noise with standard deviation 0.1. In case of Gaussian noise we use a 2 data term and in case of "salt & pepper" noise we used the 1 fidelity term. In both cases, we use TSC regularization and the curvature parameter was set to α = 10. For comparison, we also show the results of standard total variation (TV) denoising, which is obtained from our models by setting α = 0. One can clearly see that TSC regularization leads to a better preservation of image edges than TV regularization, in particular at small and elongated structures. This effect is more visible in case of "salt & pepper" noise since the 1 data term leads to a behavior similar to image inpainting, once an outlier is detected.

Proof of Theorem 1
This section and the following are devoted to the proof of Theorem 1. We first show a preliminary result which shows that admissible curves (of finite energy) in Ω × S 1 project onto curves in Ω with bounded energy as well. and assume that it is finite, which in particular yields that for a.e. s, γ x (s) = |γ x (s)|(cos γ θ (s), sin γ θ (s)) T .

Control of curves
We consider the projected curve γ x : [0, L] → Ω (possibly overlapping even if γ is simple) and wish to show that the energy E controls an energy on this curve. First we reparameterize (as usual) the rectifiable curve γ x by its length. For this we introduce the length, for s ∈ [0, L] ( (s))θ( (s))˙ (s)ds which shows that θ(l) is a BV function, or equivalently thatγ has a curvature κ =θ which is a bounded measure; in particular taking the supremum over ψ with |ψ| ≤ 1 one finds In addition, it follows, given ψ ∈ C ∞ c ( (a), (b)) with 0 ≤ a < b ≤ L that: and taking the supremum one deduces that We now show the following lemma: Lemma 5.1. Let Γ ⊂ Ω be a C 2 (oriented) curve with tangent τ Γ and curvature κ Γ , ane let and γ : [0, L] → Ω × S 1 be a rectifiable curve. Define Γ + = {x ∈ Γ : ∃s ∈ (0, L), γ x (s) = x andγ x (s) · τ Γ ≥ 0}. Then for any bounded, nonnegative, continuous function ϕ ∈ C 0 (Ω; R + ), Proof. We just need to show that, defining the measure κ as before, and the conclusion will follow from (34).
A first observation is that as τ (l) defined above is BV , it has at most a countable number of jumps and one can cover (up to the jump points) [0, (L)] with an at most countable union of intervals (a i , b i ) i∈I on which τ (l) is continuous (henceγ is a C 1 curve on (a i , b i )). Moreover, possibly dropping further a finite number of points, one may assume that |κ|(a i , b i ) < π so that γ x (a i , b i ) is a simple curve.
Assume |{l ∈ (a i , b i ) :γ(l) ∈ Γ,γ(l) · τ Γ ≥ 0}| > 0. Choose such an l ∈ (a i , b i ), with x =γ(l) ∈ Γ + , and such thatγ is differentiable in l. We assume in addition that x is a point of (H 1 -)density one inγ(a i , b i ) ∩ Γ + , that l is a Lebesgue point of κ a , the absolutely continuous part (w.r. the Lebesgue measure) of κ, and that θ(l) is approximately differentiable at l: for any η > 0, All this is true H 1 -a.e. inγ(a i , b i ) ∩ Γ + , so there is no loss of generality, see for instance [3]. Ifγ(l) = τ Γ (x), then x must be an isolated point inγ(a i , b i ) ∩ Γ + , a contradiction. Hencė γ(l) = τ Γ (x). In the same way, one has that for a.e. s near l such thatγ(s) ∈ Γ,γ(s) = τ Γ (γ(s)). For such s, one therefore has that where θ Γ is the angle between 1 0 and τΓ. Using that for such s, we obtain that for such s, (where o(1) goes to 0 as ε → 0). Hence, for ε > 0 small and η > 0, In the limit, using that l is a point of density one in the set {γ(·) ∈ Γ}, we obtain that |κ Γ (x) − κ a (l)| ≤ η and since η is arbitrary, κ Γ (x) = κ a (l). Since the intervals (a i , b i ) cover H 1 -almost all of [0, (L)], we find that this equality holds H 1 -a.e. in Γ + . The thesis of the Lemma easily follows, since

Decomposition of σ
Let u ∈ BV (Ω) with F (u) < ∞. We observe that we can decompose an admissible σ (with in particular div σ = 0) as follows: where µ is a nonnegative measure on the set of charges M which is here a shorthand notation for M(Ω × S 1 ; R 3 ): see Appendix B which discusses the results of Smirnov in [73].
Corollary 5.3. Let σ be admissible for problem (6), and assume it is decomposed as in (36).
Corollary 5.4. Let σ be admissible for problem (6), and assume it is decomposed as in (36). Let h be defined as in (3), and ϕ ∈ C 0 c (Ω) a nonnegative test function. Then Proof. This is a consequence of Lemma 5.2, however now h can take the value +∞. We simply observe that if we let, for ε > 0, then h ε is (1/ε)-Lipschitz, convex, one-homogeneous in the second variable, and continuous (and hence also (ξ, p) = ((x, θ), p) → ϕ(x)h ε (θ, p)). Hence (37) yields The results follows from the monotone convergence theorem.

Conclusion: proof of Theorem 1
We can now complete the proof of Theorem 1. We consider E a C 2 set. First, if γ(t) : [0, 1] → Ω is an oriented parameterization of (the closure of) ∂E ∩ Ω, (assuming it is connected, otherwise one needs to introduce a curve for each connected component), one can define in Ω × S 1 the curve λ(t) = (γ(t), θ(t)) where θ(t) is defined by λ(t) = |γ (t)|(cos θ(t), sin θ(t)) T . We then let Γ = λ([0, 1]) ∩ Ω × S 1 and σ = τ Γ H 1 Γ. Then, σ is admissible for (6) with u = χ E , and one finds that We need therefore to prove the reverse inequality. Let σ be admissible for (6). Thanks to [73, Theorem A] (cf Appendix B, eq. (59)), one can decompose σ as where λ are of the form λ γ = τ γ H 1 γ for rectifiable (possibly closed) curves γ ⊂ Ω × S 1 of length at most one. Equation (39) is valid for this decomposition and shows that for any ϕ ∈ C 0 c (Ω; [0, 1]), Now, for any ψ ∈ C 0 c (Ω; R 2 ), one has where τ E is a tangent vector to ∂E (oriented with E on the left-hand side and E c on the right-hand side 1 ). Let us introduce the signed distance function d E (x) = dist (x, Ω \ E) − dist (x, E) which is C 2 in a neighborhood of ∂E and is such that ∇d E = ν E (the inner normal) on ∂Ω and ∆d E = −κ E (we assume the curvature is nonnegative where the set is convex). If we consider, for ε > 0, a test function of the form Sending ε → 0 we find that Now, we can choose ϕ of the form ϕ(x)f (κ E (x)), with ϕ ∈ C 0 c (Ω; [0, 1]), assuming f is finite-valued. We obtain that Given the rectifiable curve γ, one has that thanks to (35) in Lemma 5.1. Then, (40), (42) and (43) yield and Theorem 1 is easily deduced. In case f takes the value +∞, then we first approximates f with a finite-valued function from below (replacing f with min t f (t ) + |t − t |/ε for ε > 0 small), and once the inequality is established for this function we send ε → 0, so that (9) also holds.

A Consistency of the discretization
In this appendix, we study the consistency of the discrete approximation of the problem which is used in Section 4.

A.1 Preliminary results
Let us introduce, for (s, t) ∈ R 2 ,h which is such thath (s, t) = sup as + bt.
Proof. If we introduce the (averaged) marginal σ ∈ M 1 (S 1 ; θ) defined by S 1 ψσ = 1 δ 2 x S×S 1 ψσ for all ψ ∈ C 0 (S 1 ), then one observes that This follows by a disintegration argument and using Jensen's inequality in each "slice" corresponding to a fixed value of θ. The result then follows from (44), together with Remark A.1. The last statement comes from the fact that σ x is the average of measures all contained in the cone, which is convex.
Consider now a measure σ admissible for some function u, and assume that σ is "smooth" in x: we assume for instance that it is the result of a convolution ρ ε * σ for some σ admissible (possibly extended in a larger domain), with ρ ε (x) a rotationally symmetric mollifier, as in the proof of Proposition 3.1. In this case, x → σ can be seen as a M 1 (S 1 ; R 3 )-valued smooth function.
Consider (x,θ) ∈ Ω × S 1 , δ x , δ θ small, and define in the volume V = S × [θ − δ θ /2,θ + δ θ /2) the averageσ as before and the "Raviart-Thomas" approximation of σ defined by the average fluxes through the 6 facets of V (linearly extended inside the volume, as in eq. 21). There are several ways to define this properly, at least for allθ but a countable number. In our case, one can disintegrate the measure in Ω × S 1 as σ = σ θ dµ where µ is a bounded positive measure in S 1 and for all θ, σ θ ∈ C ∞ (Ω). Then From Corollary A.5, letting e i,j,k =σ i,j,k −σ i,j,k , one has δ 2 x δ θ Proposition A.7. Let u ∈ BV (Ω), σ be admissible for u and such that Then for δ = (δ x , δ θ ) → 0 one can find a discrete field (σ 1 ) and a discrete image u δ (strongly in L 2 (Ω)) and such that for all i, j, (23) holds, for all i, j, k, (51) and (22) hold, and: whereσ is defined by (24), and whereh δ is defined in (50) (or ish in case it is Lipschitz).
To show that the discretization is consistent, we must now show a similar lower bound: namely that given any u and u δ , σ,σ which satisfy (24), (51), (23), and (52) (weakly, for instance as distributions), then one has A first obvious remark is that the field is bounded in measure, and hence, up to subsequences, converges (weakly- * ) to a measure σ. It is then easy to deduce from (51) and the convexity ofh that One can also check that div σ = 0 by passing to the limit in (22) (after a suitable integration against a smooth test function, exactly as in (55) below). Hence it is enough to show that the limiting σ is compatible with u. But this is quite obvious from (23), which one can integrate against a smooth test function, then "integrate by part" before passing to the limit. More precisely, for ϕ ∈ C ∞ c (Ω), one has (dropping the superscripts δ and denoting ϕ i,j = (1/δ 2 x ) Si,j ϕ(x)dx):  (55) as δ → 0. Eventually, we need to show a compactness property, which is that if sup δ δ 2 x δ θ i,j,kh (σ i,j,k · θ k ,σ θ i,j,k ) < +∞ (56) the discrete image u δ which is recovered from (23) (up to a constant) converges to a u(x), x ∈ Ω (in a weak sense which will be made clear). The point here is that a priori, from (56) and (1) < +∞ (57) so that the discrete total variation of u δ is a priori not well controlled. However, one can easily check that the kernel of the operator which appears in the energy (57) is two-dimensional, and made of the oscillating discrete images α, β ∈ R 2 . Hence it is possible to show that one can decompose u δ as a sum of a non-oscillating function with zero averageū δ and an oscillation v δ , and obtain a strong control on the discrete total variation ofū δ . Therefore one easily deduce that any suitably built continuous extension of u δ will converge to some u strongly in L p (Ω), for any p < 2 (as BV (Ω) is compactly embedded in such spaces), and weakly in L 2 (Ω). In addition, any control on the average of u δ and on its oscillation (which cannot be given by (57) and has to come from other terms in the energy, such as a boundary condition or a penalization: note that it is enough to control two adjacent pixels) will ensure in addition that v δ remains bounded and converges (only weakly in L p (Ω), if the control is only on the L p norm, however in this case it is obvious that the oscillating term in (58) goes to zero and v δ can only go to a constant).
To sum up, we have shown the following.
Remark A.9. In practice, we did not use the inf-convolutionsh δ (onlyh) in our discrete scheme. Also, we replaced the constraint (51) with the stronger constraintσ x i,j,k ∈ R + θ k , after having experimentally observed that there was no qualitative difference in the output. It seems the results we compute are still consistent with what is expected from the energy. B Smirnov's theorem in Ω × S 1 In this whole paper Ω ⊂ R 2 is assumed to be a Lipschitz set. In particular, locally its boundary can be represented as the subgraph {(x, y) : y < h(x)} of a Lipschitz function h. Consider a ball B where this representation holds and assume first h is C 1 , then one can extend in B a bounded Radon measure σ with div σ = 0 intoσ defined (for ψ ∈ C 0 c (B; R 2 )) Bσ · ψ := B∩Ω σ · ψ(x, y) − 1 2h (x) 0 −1 ψ(x, 2h(x) − y) .
A similar construction would allow to extend a field σ ∈ M(Ω × R 2 ; R 4 ) to M(R 4 ; R 4 ) with free divergence (either in a neighborhood or Ω × R 4 , or even everywhere). This remark allows to localize Smirnov's theorems in [73].
Consider indeed now a free divergence field σ ∈ M(Ω × S 1 ; R 3 ). It can be seen, after extension, as a field in M(R 2 × R 2 ; R 4 ) with sptσ ⊆ R 2 × S 1 .
As in Smirnov's paper [73], for l > 0 we introduce C l the set of oriented curves γ in R 4 with length l, with the topology corresponding to the weak convergence of the measures τ γ H 1 γ. Then, thanks to [73,Theorem A], σ can be decomposed as σ = C l λdµ(λ), |σ| = C l |λ|dµ(λ), for some measure µ on C l . Moreover thanks to Remark 5 in [73], µ-a.e. curve in the decomposition lies in R 2 × S 1 .
For this work, it is enough to consider l = 1. Moreover, if we restrict then all these measures to Ω × S 1 (and take for µ the corresponding marginal), we get a decomposition on curves of length less or equal to 1 (possibly entering/exiting the domain). By a slight abuse of notation we still denote C 1 such a set of curves. One finds that with now λ ∈ C 1 , the curves of length at most one in Ω × S 1 .
Remark B.1. Theorem B in [73] is a more precise statement. It shows that one can obtain a similar decomposition with now curves λ with div λ = 0 a.e.: being either finite curves entering and exiting the domain, or "elementary solenoids", which are objects of the form (in particular one should have |f (t)| = 1 a.e.), meaning that for any ϕ ∈ C 1 c (Ω × S 1 ), This expresses that either λ is defined by the closed curve f (R) (if f is periodic), or λ is a limit of curves which densify and do not loose mass in the limit. We do not need such a precise result for our construction.