extfieldsofstudy:
Computer Science
Mathematics

text:
S-MATE: Secure Coding-based Multipath Adaptive Traffic Engineering

There have been several approaches to provisioning traffic between core network nodes in Internet Service Provider (ISP) networks. Such approaches aim to minimize network delay, increase network capacity, and enhance network security services. MATE (Multipath Adaptive Traffic Engineering) protocol has been proposed for multipath adaptive traffic engineering between an ingress node (source) and an egress node (destination). Its novel idea is to avoid network congestion and attacks that might exist in edge and node disjoint paths between two core network nodes. This paper builds an adaptive, robust, and reliable traffic engineering scheme for better performance of communication network operations. This will also provision quality of service (QoS) and protection of traffic engineering to maximize network efficiency. Specifically, we present a new approach, S-MATE (secure MATE) is developed to protect the network traffic between two core nodes (routers or switches) in a cloud network. S-MATE secures against a single link attack/failure by adding redundancy in one of the operational paths between the sender and receiver. The proposed scheme can be built to secure core networks such as optical and IP networks.

I. INTRODUCTION
There have been several proposals to adapt the traffic between core network nodes in Internet Service Provider (ISP) networks [7], [12], [14]. Elwalid et al. [7] proposed an algorithm for multipath adaptive traffic engineering between an ingress node (source) and an egress node (destination) in a communication network. Their novel idea is to avoid network congestion that might exist in disjoint paths between two core network nodes. They suggested load balancing among paths based on measurement and analysis of path congestion by using Multi-Protocol Label Switching (MPLS). MPLS is an emerging tool for facilitating traffic engineering unlike explicit routing protocols that allow certain routing methodology from hop-to-hop in a network with multiple core devices. The major advantage of MATE is that it does not require scheduling, buffer management, or traffic priority in the nodes.
Network coding is a powerful tool that has been recently used to increase the throughput, capacity, and performance of wired and wireless communication networks. Information theoretic aspects of network coding have been investigated in several research papers, see for example [1], [10], [18], and the list of references therein. It offers benefits in terms of energy efficiency, additional security, and reduced delay. Network coding allows the intermediate nodes not only to forward packets using network scheduling algorithms, but also to en- Fig. 1. The network model is represented by two network nodes, an ingress node (source) and an egress node (receiver). There are k link disjoint paths between the ingress and egress nodes. code/decode them through algebraic primitive operations [1], [9], [10], [18]. For example, data loss because of failures in communication links can be detected and recovered if the sources are allowed to perform network coding operations [5], [11], [13].
MATE, which was previously proposed by one of the authors of this paper, is a traffic load balancing scheme that is suitable for S-MATE (secure MATE) as will be explained later. MATE distributes traffic among the edge disjoint paths, so as to equalize the path delays. This is achieved by using adaptive algorithms. MATE inspired other traffic engineering solutions such as TexCP [14] and the measurement-based optimal routing solution [17]. In this paper, we will design a security scheme by using network coding to protect against an entity who cannot only copy/listen to the message, but can also fabricate new messages or modify the current ones. We aim to build an adaptive, robust, reliable traffic engineering scheme for better performance and operation of communication networks. The scheme will also provide provisioning quality of service (QoS) and protection of traffic engineering.
The rest of the paper is organized as follows. In Section II, we present the network model and assumptions. In Sections III and IV, we review the MATE scheme and propose a secure MATE scheme based on network coding. Section VI provides network protection using distributed capacities and QOS, and finally Section VII concludes the paper.

II. NETWORK MODEL AND ASSUMPTIONS
The network model can be represented as follows. Assume a given network represented by a set of nodes and links. The network nodes are core nodes that transmit outgoing packets to the neighboring nodes in certain time slots. The network nodes are ingress and egress nodes that share multiple edge and node disjoint paths.
We assume that the core nodes share k edge disjoint paths, as shown in Fig. 1, for one particular pair of ingress and egress nodes. Let N = {N 1 , N 2 , ...} be the set of nodes (ingress and egress) and let L = {L 1 ℓh , L 2 ℓh , ..., L k ℓh } be the set of paths from ingress node N ℓ to an egress node N h . Every path L i ℓh carries segments of independent packets from ingress node N ℓ to egress node N h . Let P ij ℓh be the packet sent from the ingress node N ℓ in path i at time slot j to the egress N h . For simplicity, we describe our scheme for one particular pair of ingress and egress nodes. Hence, we use P ij to represent a packet in path i at time slot j.
Assume there are δ rounds (time slots) in a transmission session. For the remaining paper, rounds and time slots will be used interchangeably . Packet P ij is indexed as follows: where ID N ℓ and X ij are, respectively, the sender ID and transmitted data from N ℓ in the path L i at time slot j. There are two types of packets: plain and encoded packets. A plain packet contains the unencoded data from the ingress to egress nodes as shown in Equation (1). An encoded packet contains encoded data from different incoming packets. For example, if there are k incoming packets to the ingress node N l , then the encoded data traversed in the protection path L i lh to the egress node N h are given by where the summation denotes the binary addition. The corresponding packet becomes The following definition describes the working and protection paths between two network switches as shown in Fig. 1.
Definition 1: The working paths in a network with k connection paths carry un-encoded (plain) traffic under normal operations. The protection paths provide alternate backup paths to carry encoded traffic. A protection scheme ensures that data sent from the sources will reach the receivers in case of attack/failure incidences in the working paths.
We make the following assumptions about the transmission of the plain and encoded packets.
i) The TCP protocol will handle the transmission and packet headers in the edge disjoint paths from ingress to egress nodes. ii) The data from the ingress nodes are sent in rounds and sessions throughout the edge disjoint paths to the egress nodes. Each session is quantified by the number of rounds (time slots) n. Hence, t δ j is the transmission time at the time slot j in session δ.
iii) The attacks and failures on a path L i may be incurred by a network incident such as an eavesdropper, link replacement, and overhead. We assume that the receiver is able to detect a failure (attacked link), and our protection strategy described in S-MATE is able to recover it. iv) We assume that the ingress and egress nodes share a set of k symmetric keys. Furthermore, the plain and encoded data are encrypted by using this set of keys. That is where m i is the message encrypted by the key i . Sharing symmetric keys between two entities (two core network nodes) can be achieved by using key establishment protocols described in [15], [16]. v) In this network model, we consider only a single link failure or attack; it is thus sufficient to apply the encoding and decoding operations over a finite field with two elements, denoted as F 2 = {0, 1}. The traffic from the ingress node to the egress node in edge disjoint paths can be exposed to edge failures and network attacks. Hence, it is desirable to protect and secure this traffic. We assume that there is a set of k connection paths that need to be fully guaranteed and protected against a single edge failure from ingress to egress nodes. We assume that all connections have the same bandwidth, and each link (one hop or circuit) has the same bandwidth as the path.

III. MATE PROTOCOL
MPLS (Multipath Protocol Label Switching) is an emerging tool for facilitating traffic engineering and out-of-band control. unlike explicit routing protocols, which allow certain routing methodology from hop-to-hop in a network with multiple core devices. As shown in Fig. 2, MATE assumes that several explicit paths between an ingress node and an egress node in a cloud network have been established. This is a typical setting which exists in operational Internet Service Provider (ISP) core networks (which implement MPLS). The goal of the ingress node is to distribute traffic across the edge disjoint paths, so that the loads are balanced. One advantage of this load balancing is to equalize path delays, and to minimize traffic congestion [7], [8].
The following are the key features of the MATE algorithm. 1) The traffic is distributed at the granularity of the IP flow level. This ensures that packets from the same flow follow the same path, and hence there is no need for packet resequencing at the destination. This is easily and effectively achieved by using a hashing function on the five tuple IP address. 2) MATE is a traffic load balancing scheme, which is suitable for S-MATE, as will be explained later. MATE distributes traffic among the edge disjoint paths, so as to equalize the path delays. This is achieved by using adaptive algorithms as shown in Fig. 2 and Reference [7] 3) It is shown that distributed load balancing (for each ingress egress pair) is stable and provably convergent. MATE assumes that several network nodes exist between ingress nodes as traffic senders, and egress nodes as traffic receivers. Furthermore, the traffic can be adapted by using switching protocols such as CR-LDP [6] and RSVP-TE [4]. An ingress node is responsible to manage the traffic in the multiple paths to the egress nodes so that traffic congestion and overhead are minimized. As shown in Fig. 2, Label Switch Paths (LSPs) from the ingress node to the egress node are provisioned before the actual packet transmissions occur. Then, once the transmissions start, the ingress node will estimate the congestion that might occur in one or more of the k edge disjoint paths. As stated in [7], the congestion measure is related to one of the following factors: delay, loss rate, and bandwidth. In general, each ingress node in the network will route the incoming packets into the k disjoint paths. One of these paths will carry the encoded packets, and all other k − 1 paths will carry plain packets. Each packet has its own routing number, so that the egress node will be able to manage the order of the incoming packet, and thus achieve the decoding operations.
MATE works in two phases [7]: a monitoring phase and a load balancing phase. These two phases will monitor the traffic and balance packets among all disjoint paths. One good feature of MATE is that its load balancing algorithms equalize the derivative of delay among all edge disjoint paths from an ingress node to an egress node. Furthermore, MATE's load balancing preserves packet ordering since load balancing is done at the flow level (which is identified by a 5-tuple IP address) rather than at the packet level.

IV. S-MATE SCHEME
In this section, we provide a scheme for securing MATE, called S-MATE (Secure Multipath Adaptive Traffic Engineering). The basic idea of S-MATE can be described as shown in Equation (4). S-MATE inherits the traffic engineering components described in the previous section and in [7], [8].
Without loss of generality, assume that the network traffic between a pair of ingress and egress nodes is transmitted in rounds from ingress to egress nodes y n . . .

(4)
k edge disjoint paths, each of which carries different packets. For simplicity, we assume that the number of edge disjoint paths and the number of rounds in one transmission session are equal. There are two types of packets: i) Plain Packets: These are packets P ij sent without coding, in which the ingress node does not need to perform any coding operations. For example, in the case of packets sent without coding, the ingress node N l sends the following packet to the egress node N h : The plain data x ij are actually the encryption of the message m ij obtained by using any secure symmetric encryption algorithm [15]. That is, where key i is a symmetric key shared between N l and N h . ii) Encoded Packets: These are packets y i sent with encoded data, in which the ingress node N l sends other incoming data. In this case, the ingress node N l sends the following packet to egress node N h : The encoded packet will be used in case any of the working paths is compromised. The egress node will be able to detect the compromised data, and can recover them by using the data sent in the protection path. Lemma 2: The S-MATE scheme is optimal against a single link attack.
What we mean by optimal is that the encoding and decoding operations are achieved over the binary field with the least computational overhead. That is, one cannot find a better scheme than this proposed encoding scheme in terms of encoding operations. Indeed, one single protection path is used in case of a single attack path or failure. The transmission is done in rounds (time slots), and hence linear combinations of data have to be from the same round time. This can be achieved by using the time slot that is included in each packet sent by the ingress node.
Lemma 3: The network capacity between the ingress node and the egress node is given by k − 1 in the case of one single attack path. node node y 1 x 21 x 31 x 41 x 51 x 11 y 2 x 32 x 42 x 52 x 12 x 22 y 3 x 43 x 53 x 13 x 23 x 33 y 4 x 54 x 14 x 24 x 34 x 44 y 5 Fig. 3. Working and protection edge disjoint paths between two core nodes. The protection path carries encoded packets from all other working paths between ingress and egress nodes.
Encoding Process: There are several scenarios where the encoding operations can be achieved. The encoding and decoding operations will depend mainly on the network topology; i.e., on how the senders and receivers are distributed in the network.
• The encoding operation is performed at only one ingress node N l . In this case, N l will prepare and send the encoded data over L i lh to the receiver N h . • We assume that k packets will be sent in every transmission session from the ingress node. Also, if the number of incoming packets is greater than k, then a mod function is used to moderate the outgoing traffic in k different packets. Every packet will be sent in a different path. • Incoming packets with large sizes will be divided into small chunks, each with an equal size. Decoding Process: The decoding process is performed in a similar way as explained in the previous works [2], [3].
We assume that the ingress node N l assigns the paths that will carry plain data as shown in Fig. 3. In addition, N l will encode the data from all incoming traffics and send them in one path. This will be used to protect any single link attacks/failures. The objective is to withhold rerouting signals or transmitted packets due to link attacks. However, we provide strategies that use network coding and reduced capacity at the ingress nodes. We assume that the source nodes (ingress) are able to perform encoding operations and the receiver nodes (egress) are able to perform decoding operations.
One of S-MATE's objectives is to minimize the delay of the transmitted packets. So, the packets from one IP address will be received in order on one path. The key features of S-MATE can be described as follows: • The traffic from the ingress node to the egress node is secured against eavesdroppers and intruders. • No extra paths in addition to the existing network edge disjoint paths are needed to secure the network traffic. • It can be implemented without adding new hardware or network components. The following example illustrates the plain and encoded data transmitted from five senders to five receivers. Example 4: Let N l and N h be two core network nodes (sender and receiver) in a cloud network. Equation (7) explains the plain and encoded data sent in five consecutive time slots from the sender to the receiver. In the first time slot, the first connection carries encoded data, and all other connections carry plain data. Furthermore, the encoded data are distributed among all connections in the time slots 2, 3, 4 and 5.
The encoded data y j , for 1 ≤ j ≤ 5, are sent as We notice that every message has its own time slot. Hence, the protection data are distributed among all paths for fairness.

V. A STRATEGY AGAINST TWO ATTACKED PATHS
In this section, we propose a strategy against two attacked paths (links), securing MATE against two-path attacks. The strategy is achieved using network coding and dedicated paths. Assume we have n connections carrying data from an ingress node to an egress node. All connections represent disjoint paths.
We will provide two backup paths to secure against any two disjoint paths, which might experience any sort of attacks. These two protection paths can be chosen using network provisioning. The protection paths are fixed for all rounds per session from the ingress node to the egress node, but they may vary among sessions. For example, the ingress node N l transmits a message x iℓ to the egress node N h through path L i ℓh at time t ℓ δ in round time ℓ in session δ. This process is explained in Equation (9)  x n1 x n2 x n3 . . . x nn . . . 5 All y ℓ j 's are defined as: The coefficients a ℓ i and b ℓ i are chosen over a finite field F q with q > n − 2, see [3] for more details. One way to choose these coefficients is by using the follow two vectors.
Therefore, the coded data is In the case of two failures, the receivers will be able to solve two linearly independent equations in two unknown variables. For instance, assume the two failures occur in paths number two and four. Then the receivers will be able to construct two equations with coefficients Therefore, we have One can multiply the first equation by α and subtract the two equations to obtain value of x 4ℓ . We notice that the encoded data symbols y jℓ and y kℓ are fixed per one session but it is varied for other sessions. This means that the path L j lh is dedicated to send all encoded data y j1 , y j2 , . . . , y jn .
Lemma 5: The network capacity of the protection strategy against two-path attacks is given by n − 2.
There are three different scenarios for two-path attacks, which can be described as follows: i) If the two-path attacks occur in the backup protection paths L j lh and L k lh , then no recovery operations are required at the egress node. ii) If the two-path attacks occur in one backup protection path say L j lh and one working path L i lh , then recovery operations are required. iii) If the two-path attacks occur in two working paths, then in this case the two protection paths are used to recover the lost data. The idea of recovery in this case is to build a system of two linearly independent equations with two unknown variables.

VI. NETWORK PROTECTION USING DISTRIBUTED CAPACITIES AND QOS
In this section, we develop a network protection strategy in which some connection paths have high priorities (less bandwidth and high demand). Let k be the set of available connections (disjoint paths from ingress to egress nodes). Let m be the set of rounds in every cycle. We assume that all connection paths might not have the same priority demands and working capacities. Connections that carry applications with multimedia traffic have higher priorities than those of applications carrying data traffic.. Therefore, it is required to design network protection strategies based on the traffic and sender priorities. Consider that available working connections k may use their bandwidth assignments in asymmetric ways. Some connections are less demanding in terms of bandwidth requirements than other connections that require full capacity frequently. Therefore, connections with less demand can transmit more protection packets, while other connections demand more bandwidth, and can therefore transmit fewer protection packets throughout the transmission rounds. Let m be the number of rounds and t δ i be the time of transmission in a cycle δ at round i. For a particular cycle i, let t be the number of protection paths against t link failures or attacks that might affect the working paths. We will design a network protection strategy against t arbitrary link failures as follows. Let the source s j send d i data packets and p i protection packets such that d j + p j = m. Put differently: In general, we do not assume that d i = d j and p i = p j . The encoded data y iℓ are given by We assume that the maximum number of attacks/failures that might occur in a particular cycle is t. Hence, the number of protection paths (paths that carry encoded data) is t. The selection of the working and protection paths in every round is performed by using a demand-based priority function at the senders's side. It will also depend on the traffic type and service provided on these protection and working connections. See Fig. 4 for ingress and egress nodes with five disjoint connections.
In Equation (17), every connection i is used to carry d i unencoded data x i1 , x i2 , . . . , x idi (working paths) and p i Fig. 4. Working and protection edge disjoint paths between two core nodes (ingress and egress nodes). Every path L i carries encoded and plain packets depending on the traffic priority p i and time t i . encoded data y i1 , y i2 , . . . , y ipi (protection paths) such that d i + p i = m.
Lemma 6: Let t be the number of connection paths carrying encoded data in every round. Then, the network capacity C N is given by Proof: The proof is straightforward from the fact that t protection paths exist in every round, and hence k − t working paths are available throughout all m rounds.

VII. CONCLUSION
In this paper, we have proposed the S-MATE scheme (secure multipath adaptive traffic engineering) for operational networks. We have used network coding of transmitted packets to protect the traffic between two core nodes (routers, switches, etc.) that could exist in a cloud network. Our assumption is based on the fact that core network nodes share multiple edge disjoint paths from the sender to the receiver. S-MATE can secure network traffic against single link attacks/failures by adding redundancy in one of the operational paths. Furthermore, the proposed scheme can be built to secure operational networks including optical and multipath adaptive networks. In addition, it can provide security services at the IP and data link layers.


extfieldsofstudy:
Mathematics
Computer Science
Physics

text:
Spatially regularized compressed sensing of diffusion MRI data

The present paper introduces a method for substantial reduction of the number of diffusion encoding gradients required for reliable reconstruction of HARDI signals. The method exploits the theory of compressed sensing (CS), which establishes conditions on which a signal of interest can be recovered from its under-sampled measurements, provided that the signal admits a sparse representation in the domain of a linear transform. In the case at hand, the latter is defined to be spherical ridgelet transformation, which excels in sparsifying HARDI signals. What makes the resulting reconstruction procedure even more accurate is a combination of the sparsity constraints in the diffusion domain with additional constraints imposed on the estimated diffusion field in the spatial domain. Accordingly, the present paper describes a novel way to combine the diffusion- and spatial-domain constraints to achieve a maximal reduction in the number of diffusion measurements, while sacrificing little in terms of reconstruction accuracy. Finally, details are provided on a particularly efficient numerical scheme which can be used to solve the aforementioned reconstruction problem by means of standard and readily available estimation tools. The paper is concluded with experimental results which support the practical value of the proposed reconstruction methodology.

signals. What makes the resulting reconstruction procedure even more accurate is a combination of the sparsity constraints in the diffusion domain with additional constraints imposed on the estimated diffusion field in the spatial domain. Accordingly, the present paper describes a novel way to combine the diffusion-and spatial-domain constraints to achieve a maximal reduction in the number of diffusion measurements, while sacrificing little in terms of reconstruction accuracy. Finally, details are provided on a particularly efficient numerical scheme which can be used to solve the aforementioned reconstruction problem by means of standard and readily available estimation tools. The paper is concluded with experimental results which support the practical value of the proposed reconstruction methodology.

Introduction
The human brain consists of about 10 11 nerve cells that can be subdivided into about 1000 different cell types, a complexity that far exceeds that of other organs of the body. A further complexity is evident in the way in which the component cells of the brain interconnect and function [1]. In contrast to other types of the cells, each neuron communicates with many target cells by means of its distinctive protoplasmatic protrusion, called an axon. Axons with similar destinations, in turn, tend to form bundles -known as neural fibre tractswhich play a pivotal role in the determination of brain connectivity. Through reconstructing the pattern of connectivity of the neural tracts in both healthy and diseased subjects, it is therefore possible to obtain an abundance of valuable diagnostic information that could be used for early diagnostics of brain-related disorders, for assessing the damage caused to the brain by stroke, tumours or injuries, as well as for planning and monitoring of neurosurgeries [2].
Central to MRI is the notion of contrast, which is typically defined by the biochemical composition of interrogated tissue as well as by the morphology of its associated parenchyma. Prevalent in MRI practice are the contrasts determined by the T 1 /T 2 relaxation times and proton density (PD). Despite their exceptional importance to clinical diagnosis, none of the above contrast mechanisms has demonstrated effectiveness in delineating the morphological structure of the white matter. It is only with the advent of diffusion MRI (dMRI) that scientists have been able to perform quantitative measurements of the diffusivity of white matter, based on which its structural delineation has become possible [2][3][4][5][6][7][8]. Although over the last two decades dMRI has developed into an established technique with a great impact on health care and neurosciences, like any other MRI technique it remains subject to artifacts and pitfalls [3]. While many of such artifacts can be overcome by using advanced hardware designs and/or more sophisticated imaging protocols [2,Ch.2], [9], one particularly critical limitation of dMRI stems from the physics of the acquisition of diffusion MR images, and therefore is impossible to resolve by operational means. Specifically, since collecting the dMRI data requires a repetitive acquisition of MR responses from the same volume for a number of diffusion-encoding gradi-ents, it is the relatively long acquisition times that greatly impair the practical value of this important imaging modality. Particularly, longer acquisition times entail a higher probability for the patient to exercise involuntary motion (typically caused by fatigue and/or stress related tremors, swallowing, uncontrollable sighing or coughing), which severely affects the quality of dMRI data. Moreover, since during the whole duration of the scan the patients are required to remain motionless, it is currently deemed ineffective to apply dMRI-based diagnosis in paediatrics as well as to patients with dementia or post-traumatic injuries, where non-compliance is typical. The problem of long acquisition times also hampers the application of dMRI for intra-operational imaging, where it could be an irreplaceable tool to use for neurosurgical planning and decision-making support [10][11][12]. Lastly, relatively long scanning times required by dMRI aggravate the problem of accessibility to MRI equipment. All the above arguments suggest that the practical value of dMRI could be improved by shortening the scanning times required for acquisition of dMRI data. A particular method to achieve such an improvement is detailed in this paper.
In this work, we adopt a general diffusion model in which each voxel within a region of interest (ROI) is allowed to support more than one fibre tract. In this case, under some general assumptions (see, e.g., [13,Sect. 3.1] for more details), the diffusion signal s(u | r) originating from a voxel with spatial coordinate r ∈ R 3 containing M (r) fibres can be modelled as [13,14] s(u | r) = s 0 (r) where u denotes the spherical coordinate, i.e.
u ∈ S 2 := v ∈ R 3 | v 2 = 1 , and α i (r) > 0 are positive weights obeying M (r) i=1 α i (r) = 1. In (1), s 0 denotes the diffusion signal obtained in the absence of diffusion encoding (i.e. the socalled "b0 image"), b is defined as a function of the shape and amplitude of diffusion-encoding gradients [15,Eq. 3.18], and {D i (r)} M (r) i=1 are 3 × 3 diffusion tensors associated with the M (r) neural fibres passing through the coordinate r. In practical settings, the spherical coordinate u is sampled at K distinct points {u k } K k=1 over the unit sphere. In this case, for each u k , MR measurements are acquired in the form of a diffusion-encoded image s k (r) := s(u k | r). As a result, a typical dMRI data set consists of a collection of such diffusion-encoded images {s k (r)} K k=1 , whose size K determines the accuracy with which the directions of local diffusion flows can be estimated.
Provided unlimited scanning time, one could measure the diffusion in thousands of orientations, making it possible to identify the directions of dominant diffusion with very high precision. For the reasons explained earlier, however, scanning times are always limited, which necessitates restriction of K to a reasonably small value. This brings us to the central question addressed in the present paper: what is a sufficient number K of diffusion-encoding directions to use? It turns out that, for some realizations of dMRI, the above question can be answered in a rigorous manner. In particular, in diffusion tensor imaging (DTI) [3-5, 7, 8, 16, 17], the reconstruction is carried out under the assumption that each voxel can support only one diffusion flow as most, which corresponds to setting M (r) = 1 for all r in (1). Accordingly, a minimum of K = 7 diffusionencoded images are theoretically sufficient to measure S 0 (r) and recover the six non-repetitive components of the symmetric tensors D(r) by means of leastsquare fitting. (In practice, however, a larger number of gradient directions is employed to reduce the estimation variance, with a typical K being between 25 and 30 [18].) Unfortunately, the accuracy of DTI is known to deteriorate dramatically at the sites where the neural fibres (or bundles thereof) cross, touch upon each other, or diverge [13,[19][20][21][22][23][24][25][26][27].
The fibre crossing problem in DTI has prompted efforts to develop dMRI methodologies which are capable of detecting multiple diffusion flows (or, equivalently, neural fibre tracts) within a given voxel. One of such techniques is High Angular Resolution Diffusion Imaging (HARDI) [14,20,[22][23][24][26][27][28], which is capable of capturing multi-modal diffusion patterns by sampling the spherical shell at a much greater number of points (usually between 60 and 100) as compared to the case of DTI. Increasing K makes it possible to describe the diffusion measurements using much more accurate models. Among these are parametric models [29][30][31][32][33][34] which allow HARDI signals to be expressed in terms of a relatively small number of prototype signal forms. Unfortunately, fitting a parametric model often entails minimization of non-convex functionals, which is a noise-sensitive and computationally intensive task, prone to the problem of local minima. The need to predetermine the optimal number of fitting terms is known to be another disadvantage of using the models of the above type.
The problems associated with parametric modelling of HARDI signals can be overcome by using non-parametric models, in which case the signals are recovered by projecting the observed data onto properly defined functional subspaces. In particular, the applicability of spherical Fourier analysis [35] to dMRI has been demonstrated in [25][26][27][28], where HARDI signals are approximated by truncated series of spherical harmonics (SH). Despite its stability and computational efficiency, however, the SH-parameterization involves a relatively large number of SHs, which suggests that the SHs cannot be an adequate basis for sparse representation of HARDI signals. The main reason for this is rooted in the fact that the energies of elementary signals (1) are concentrated alongside the great circles of S 2 , whereas the energy of SHs is spread all over S 2 , and, as a result, a relatively large number of SHs are needed to effectively "encode" each of d i . The inability of the basis of SHs to sparsely represent diffusion signals has led to the proposal of spherical ridgelets in [36,37], where it was shown that it only takes 6 to 8 spherical ridgelets on average to represent the HARDI signals with a precision exceeding the precision of their representation with 45 SHs.
The present work takes the ideas of [36,37] one step further and shows that the availability of a sparsifying basis for HARDI signals can be used to reduce the number of diffusion gradients required for data acquisition. In particular, we suggest to use the theory of compressed sensing (CS) [38][39][40][41][42][43][44] to recover the HARDI signals using the number of spherical samples K in a range of values typical for DTI (i.e. K ∈ [16,24]), thus allowing a multi-fibre analysis of dMRI data to be performed at the "acquisition cost" of a standard DTI.
It is worthwhile noting that the ideas of CS have already paved their way into the field of diffusion imaging [45][46][47][48][49]. In this regard, conceptually the closest to the proposed approach is the method reported in [46]. In spite of this similarity, however, there are two principal distinctions which make the present method a more powerful alternative. In particular, the basis functions used in [46] are limited to represent an average diffusivity and anisotropy of the white matter, thereby neglecting both intra-and inter-voxel variability of tensors D i (r) in (1). The ridgelet representation, on the other hand, is a multiresolution technique, which possesses an intrinsic ability to deal with a continuum of different diffusion scales. Second, the approach in [46] is applied in a "voxel-by-voxel" manner and it therefore does not take into consideration the spatial regularity of diffusion field. The present paper, on the other hand, proposes a novel formulation of the problem of CS-based reconstruction of diffusion signals, in which the sparsity constraints enforced in the diffusion domain are augmented by regularity constraints enforced in the spatial domain. The resulting reconstruction problem has the format of a convex minimization problem, which is solved using a specially adapted version of the split Bregman algorithm [50][51][52]. As will be shown below, the proposed algorithm results in a particularly advantageous computational structure which allows the solution to be computed via a sequence of simple and easily parallelizable steps.
The rest of the paper is organized as follows. Section II provides additional comments on the input-output structure of the proposed algorithm. The construction of spherical ridgelets is briefly outlined in Section III, whereas Section IV gives a formal description of the proposed reconstruction methodology. Some principal details on the numerical implementations of the proposed algorithm are summarized in Section V, with the results of our experimental studies reported in Section VI. Section VII finalizes the paper with a discussion and conclusions.

Problem Statement
In the centre of our considerations is the diffusion signal s(u | r) which, when normalized by its related b0-image s 0 (r), quantifies the attenuation of MR readout caused by the diffusion of water molecules in the direction u ∈ S 2 through the spatial position r ∈ R 3 . In practical settings, both u and r are discretized. Specifically, restricting u to a discrete set of orientations {u k } K k=1 prescribes the acquisition of diffusion data in the form of K diffusion-encoded images {s k (r)} K k=1 , with each s k (r) corresponding to a given u k . In this case, for a fixed r 0 , the vector [s 1 (r 0 ), s 2 (r 0 ), . . . , s K (r 0 )] T ∈ R K represents a discretization of s(u | r 0 ). Note that such a discretization follows a linear measurement model, since each sample s k (r 0 ) can be represented as an inner product of s(u | r 0 ) with a sampling function. In particular, let {ϕ k (u)} K k=1 be a Dirac basis of sampling functions defined by where δ denotes the Dirac delta function and the dot stands for the Euclidean dot product. Then, formally, with dη being the standard rotation invariant measure on S 2 .
Next, given a collection of M spherical ridgelets {ψ m (u)} M m=1 (defined below), the signal s(u | r) is assumed to be expandable as with c(r) ∈ R M being a vector of spherical ridgelet coefficients which depend on the spatial coordinate r. It is important to note that the set of spherical ridgelets is allowed to be overcomplete, implying dim Span{ψ m (u)} M m=1 < M . A practical consequence of this fact is that the definition of coefficients c(r) in (5) is, in general, not unique. This non-uniqueness is further aggravated by the fact that c(r) will have to be recovered from an under-sampled set of diffusion measurements, in which case K M . Overcoming such a severe underdetermination in the problem of estimating the ridgelet coefficients c(r) will be possible based on the fundamental premise of the theory of CS, which states that an accurate estimation of c(r) is possible if the latter is sufficiently sparse and if the sampling and representation bases are sufficiently decorrelated. While the sparsity of c(r) is rooted in the very design of spherical ridgelets [37], the incoherency between the Dirac sampling functions (3) and spherical ridgelets stems from the fact that the former have an infinitely small support, whereas the latter are "smeared" all over the unit sphere. The above properties of the basis of spherical ridgelets yield conditions for an effective application of CS, in which case one can obtain a faithful reconstruction of diffusion signals using as few as K = 16 diffusion-encoding gradients.
The proposed algorithm produces an estimate of the ridgelet representation coefficients c(r) in (5). Once available, the coefficients provide an access to the analytical definition of diffusion signals by virtue of (5). This can be used in a number of ways. One possibility could be to use the ridgelet coefficients to compute their associated orientation distribution functions [14], based on which a multi-fibre tractography analysis can be done [53,54]. Alternatively, (5) can be used to evaluate the diffusion signals over an arbitrarily fine grid of orientations. Subsequently, such refined "measurements" could be fitted using a different representation model, whose application to the original data would not have been possible without causing severe underestimation errors. Deconvolving the refined data to estimate the underlying fibre orientation functions [55][56][57][58] would be another important option to follow. In this paper, we refrain from questioning which of the above possibilities is more advantageous over the others. Our sole objective here is to specify a signal processing algorithm which can be used to recover HARDI signals, while using the number of diffusion-encoded images typical for a standard DTI.
Finally, it should be noted that the primarily purpose of the proposed methodology is to improve the value of HARDI in terms of its time efficiency. Since the improvement is achieved through merely decreasing the number of diffusion-encoding gradients, the proposed method by no means abrogates the use of fast imaging protocols [59][60][61] to further accelerate the data acquisition. Furthermore, an additional speed-up can be achieved via applying CS to reconstruct the diffusion encoded images s k (r) from their subcritical samples in the spectral domain [62][63][64][65][66]. Generally speaking, we believe this is a combination of such software and hardware technologies which will eventually lead to substantantial improvements in the practical value of HARDI-based diagnosing. In this paper, however, we confine our contribution to showing one particular way of attaining this important objective.

Spherical ridgelets
It is the property of spherical ridgelets to provide sparse representation of diffusion signals described by (1) which makes them an unparalleled tool for CSreconstruction of HARDI data. To avoid repetitions, in what follows, we present only the most principal points of ridgelets design, while their detailed description can be found in [37].
Spherical ridgelets are constructed using the fundamental principles of wavelet theory [67,68]. Specifically, let x ∈ R + and ρ ∈ (0, 1) be a positive scaling parameter. Further, let κ(x) = exp{−ρ x (x + 1)} be a Gaussian function, which we subject to a series of dyadic scalings to result in where j ∈ N := {0, 1, 2, . . .}. Subsequently, the Gaussian-Weierstrass scaling function χ j,v : S 2 → R at resolution j ∈ N and orientation v ∈ S 2 can be defined as given by [69,70] where P n denotes the Legendre polynomial of order n. It is worth noting that the L 2 energy of χ j,v is concentrated around the spherical point v, with this concentration becoming more and more localized when j approaches infinity. The spherical ridgelets are designed with the help of the Funk-Radon transform which, for an arbitrary continuous function f : S 2 → R, is defined as with σ(v) denoting the great circle perpendicular to direction v, i.e. σ(v) := u ∈ S 2 | u · v = 0 . Subsequently, following [37], the semi-discrete frame U of spherical ridgelets can be defined as where the spherical ridgelet functions ψ j,v are obtained from χ j,v according to with χ −1,v (u) ≡ 0. Using (7), the ridgelets (10) can be redefined in a closed form as (see [35] for details) where κ −1 (n) = 0, ∀n and The set U in (9) is infinite-dimensional, and hence is not suitable for practical computations. To define a discrete counterpart of U, one has first to restrict the values of the resolution index j to a finite set {−1, 0, 1, . . . , J}, where J defines the highest level of "detectable" signal details. Additionally, the set of all possible orientations v ∈ S 2 of spherical ridgelets needs to be discretized as well. To find a proper discretization scheme, we first note that the construction in (11) suggests that the bandwidth of the spherical ridgelets (and therefore the dimensionality of the functional space they belong to) increases proportionally to 2 j . Since the space of spherical harmonics of degree n has a dimension of (n + 1) 2 , it seems to be reasonable to define the number of ridgelet orientations at resolution j to be equal to M j = (2 j+1 m 0 + 1) 2 , with m 0 being the smallest spherical order resulting in κ 0 (m 0 ) ≤ for some predefined 0 < 1 (e.g. = 10 −6 ). Consequently, for each j, a total of M j orientations {v i j } Mj i=1 are chosen so that a discrete counterpart of U can now be defined as where the subscript d stands for "discrete". It should be noted that, although the set U d is composed of continuously defined functions, its dimension is finite, since U d consists of a total of M = J j=−1 (2 j+1 m 0 + 1) 2 spherical ridgelets. To slightly simplify our notation, in what follows, the spherical ridgelets in U d will be indexed as ψ m (u), with m = 1, 2, . . . , M being a combined index accounting for both different resolutions and orientations.
Given a sampling set of K diffusion-encoding orientations {u k } K k=1 , one can use (11) to compute the values of the spherical ridgelets in U d over the sampling set 1 . The resulting values can be stored into a K × M matrix A defined as Then, for a given vector s(r) : T ∈ R K of the measured values of a diffusion signal s(u | r) at the spatial location r, the model (5) asserts the existence of representation coefficients c(r) ∈ R M such that where e(r) accounts for both model and measurement noises. Clearly, the nonnegligibility of e(r) along with the fact that K M makes the problem of recovering the representation coefficients c(r) from s(r) a very challenging inverse problem, our solution to which is presented next.

Proposed reconstruction framework
Let Ω represent the volume within which diffusion measurements are acquired. Also known as a region of interest, Ω is assumed to be a bounded rectangular subdomain of R 3 , i.e. Ω : Let Ω d be a discrete subset of Ω, which represents the spatial locations at which the diffusion signal is measured. Specifically, Ω d is assumed to be a uniform lattice which can be formally defined as where 0 ≤ i 1 < N x , 0 ≤ i 2 < N y , and 0 ≤ i 3 < N z are sampling indices in the direction of x, y and z coordinates, respectively. Let K be the number of diffusion-encoding gradients used for HARDI data acquisition, and let the corresponding gradient orientations be denoted by {u k } K k=1 , where u k ∈ S 2 . For each of these values of u k , MRI measurements result in its associated diffusion-encoded image s k , which can be formally viewed as a mapping from Ω d to R. For the sake of convenience, each s k can be stored and manipulated as an N x × N y × N z array of real numbers, namely s k ∈ R Nx×Ny×Nz . Alternatively, at a given coordinate r ∈ Ω d , one can combine the values s 1 (r), s 2 (r), . . . , s K (r) into a column vector s(r) := [s 1 (r), s 2 (r), . . . , s K (r)] T ∈ R K (as it was already done in (15)). This vector can then be regarded as a vector of discrete measurements of an associated HARDI signal s(u | r) measured at orientations {u k } K k=1 . It is worth noting that, according to the above notations, the value s k (r) admits a twofold interpretation, viz. either as the k th coordinate of vector s(r) or as the value of image s k at spatial position r.
When combined together, the continuum of vectors s(r) can be regarded as a discrete vector field s : Ω → R K , in which case s(r) has a natural interpretation of the value of s corresponding to position r. The vector space V of such vector fields can be endowed with the standard inner product being the standard inner product between the scalar-valued images s 1 k and s 2 k . Accordingly, congruent to the definition in (17), where · 2 and · F denote the Euclidean vector and Frobenius matrix norms, respectively.
Another norm on V that we shall make use of is the total variation (TV) semi-norm which is defined as follows. First, let us define the total variation of the k th component s k of the field s in a standard manner as where Consequently, the TV norm of the discrete vector field s can be defined in terms of the TV-norms of its K components as Thus, for example, α = 2 was used in the TV-denoising method reported in [71].
In this paper, we use α = 1. Now, let {ψ m } M m=1 be a set of spherical ridgelets defined by (13), which is assumed to be rich enough so that each HARDI signal can be expressed according to (5). Analogously to the discrete measurements s(r), the representation coefficients c(r) corresponding to different voxels r can be aggregated into a vector field c ∈ U, where U : Ω d → R M (with c(r) being the value of c observed at r). Although it is possible to endow the vector space U with both the 2and TV-norms by analogy with (18) and (20), it will be particularly useful to consider the 1 -norm of c which can be defined as Using the definitions of the vector fields V and U as well as the definition of A in (14), a connection between V and U is established by means of a linear map A : U → V that is defined as given by Consequently, using A, one can define the HARDI data formation model as where e ∈ V is supposed to account for both measurement noise and modelling errors, and it is assumed to have a relatively small 2 -norm e V,2 ≤ . The model (23) suggests a reduction of the problem of estimation of HARDI signals to the problem of estimation of their corresponding representation coefficients c from the discrete and noisy measurements s. Furthermore, as our main intension has been to recover the coefficients c using as few diffusion-encoding gradients as possible (implying K M ), there is an infinite number of solutions which would fit the constraint A{c} − s V,2 ≤ . However, if it is known a priori that, for each r ∈ Ω d , the vector of representation coefficients c(r) is sparse, then a useful estimate of c can be obtained as a solution to the following convex optimization problem [38][39][40][41][42][43][44] It should be noted that the optimization problem (24) is equivalent to solving and therefore, under the assumption of spatially homogeneous noise e, the problem (24) is separable in the spatial coordinate r. This means that an optimal field c can be recovered by solving for its components min c(r) independently at each r ∈ Ω d . While computationally attractive, the above solution is suboptimal, since it completely disregards the dependencies which are likely to exist between spatially adjacent HARDI signals. A possible way to take such dependencies into consideration is to require the noise-free version of the measured signal s to possess a minimal TV norm among all possible candidate solutions [72]. This requirement can be translated into the following minimization problem where the role of γ > 0 is to balance the relative influence of the sparse and TV terms in the above cost function. The optimization problem (30) can be rewritten in its equivalent Lagrangian form for some optimal values of λ > 0 and µ > 0 [73]. Below, we are going to specify a particular, computationally efficient method for solving (31). In this connection, it is instructive to outline the following two instances of (31).

Sparse-only reconstruction
When µ = 0, the functional in (31) becomes separable in the spatial variable r in the sense that, in such a case, an optimal c can be recovered by solving min c(r) for each c(r) independently. Note that (32) can be considered to be a Lagrangian form of the optimization problem (28). There exist a broad spectrum of methods which could be used for solution of (32). Some particularly attractive algorithms seem to be those exploiting the principle of iterative shrinkage (aka iterated thresholding) [74][75][76][77]. While the non-differentiability of 1 -norm in (32) rules out the applicability of gradient-based optimization tools, iterative shrinkage methods are capable of finding a solution of (32) through iterative application of a first-order, fixed-point update rule. Specifically, many algorithms of this kind find an optimal solution as a stationary point of the sequence of estimates produced by where S τ {t} = sign(t)(|t| − τ ) + denotes the operator of soft thresholding and ν is chosen to obey ν > AA T . In the present paper, a modification of the iterative update in (33), known as the fast iterative shrinkage thresholding algorithm (FISTA) [78], was employed due to its considerably faster convergence as compared with many alternative "accelerated" methods. It should be emphasized that, while being suboptimal from the viewpoint of spatial-domain regularity, the solution of (32) through iterative shrinkage is advantageous in two important practical ways. First, it suggests considerable storage reduction, since the thresholding operator in (33) sets to zero the representation coefficients with amplitudes less or equal to λ/ν in absolute value. It makes it possible to use sparse data formats to store and manipulate the representation coefficients. Second, the fact that the estimation of c(r) is performed at each voxel independently suggests a natural way to speed up the overall estimation process though parallel computing on a multicore system.

TV-only reconstruction
When λ = 0, solving the optimization problem (31) is equivalent to simultaneously solving K optimization problems of the form where k = 0, 1, . . . , K − 1 and [A{c}] k denotes the k-th component of the vector field A{c} ∈ V. Let [A{c}] k be denoted by u k , i.e. u k := [A{c}] k . Then, reformulated with respect to u k , the problem (34) can be rewritten as in which case it can be recognized as the problem of TV-denoising of the diffusion-encoded image s k [72]. It is important to note that, in contrast to (34), the problem (35) can be solved for each k independently, in which case we say that the estimation becomes separable in the diffusion direction.
The current arsenal of methods which can be used for solving (35) is broad. Originated from the work of Rudin et al [72], TV-denoising methods currently include gradient-based optimization methods [79], first-order methods [80], iterative shrinkage methods [81,82], and Bregman-type iterative algorithms [50,51,83]. The first-order methods appear to be a particularly attractive option when relatively large data arrays need to be processed, which is the case relevant to dMRI. Consequently, in the present paper, we employ the algorithm of [84] for the simplicity and elegancy of its implementation as well as for its outstanding convergence properties.

Solution using Split Bregman Algorithm
Directly solving the original problem (31) is difficult because of the compound nature of the regularization it involves. The split Bregman approach [50,51] allows reducing (31) to a simpler form through introduction on an auxiliary variable u ∈ V, which can be viewed as a noise-free version of the data field s. Particularly, using u (31) can be redefined as min c,u s.t. A{c} − u 2 V,2 = 0.
Then, starting from an arbitrary b 0 ∈ V, the Bregman algorithm [85] finds optimal c and u through the following iterations u t+1 , c t+1 = arg min c,u for some γ > 0 2 . The functional in (37) is supposed to be minimized over two variables, i.e. u and c. However, due to the way the 1 and TV components of this functional have been split, the minimization can now be performed efficiently by iteratively minimizing with respect to u and c separately. The resulting iteration steps are Step 1: c t+1 = arg min Step 2: u t+1 = arg min Note that the functional at Step 2 contains two quadratic terms which can be combined together to result in Step 2: u t+1 = arg min To cause a substantial reduction in the value of the cost functional in (37), Step 1 and Step 2 should be applied recursively for a predefined number of times before the Bregman parameter b t is updated according to (37). It was argued in [51], however, that the extra precision gained through such a repetitive application of Step 1 and Step 2 is likely to be "wasted" when b t is updated. Consequently, it was suggested in [51] to perform these steps only once per iteration cycle. It is interesting to note that, in this case, the split Bregman algorithm transforms into the alternating directions method of multipliers (ADMM) [52], whose convergence is guaranteed by the Eckstein-Bertsekas theorem [86] (see also Theorem 3.1 in [52]).
The final algorithm is summarized below. Lines 3-4 of Algorithm 1 correspond to Step 1 in (38), while lines 5-6 correspond to Step 2. An even more important fact to notice is that the optimization problem in line 4 is separable in the spatial coordinate r. This optimization, therefore, can be performed at each voxel independently as discussed in Section 4.1. Moreover, the optimization problem in line 6 is separable in the diffusion coordinate k, and hence it amounts to applying TV-denoising to each of the K components of u independently as discussed in Section 4.2.

Technical details of the experimental study
Both the choice of diffusion-encoding directions {u k } K k=1 and of the orientations of spherical ridgelets require the use of a sampling scheme. In dMRI, one of the standard methods used to distribute a given number of spherical points over S 2 in a quasi-uniform manner is by means of the electrostatic repulsion algorithm (also known as the Thomson problem). However, since the diffusion signals are symmetric (implying s(u | r) = s(−u | r)), it is a unit hemisphere, not the entire S 2 , which actually needs to be discretized. In view of the absence of a formulation of the Thomson problem for hemispherical domains, a common practice is to run the standard procedure for twice as many points as needed, followed by keeping only a half of the resulting configuration. However, as the retained points are not explicitly constrained to lie on a hemisphere, they may include nearly antipodal pairs which are likely to introduce undesirable dependencies between the diffusion measurements as well as between the basis functions. This limitation can be overcome through adapting a different sampling strategy. Particularly, in this paper, both the diffusion-encoding directions and the orientations of spherical ridgelets have been defined by using the method of generalized spiral points [87,88], in which the sampling points are arranged along a spiral in such a way that the distance between the points along the spiral is approximately equal to the distance between its coils. This method is easily adaptable for sampling of the "northern" hemisphere (i.e. u ∈ S 2 | u · [0, 0, 1] T ≥ 0 ), providing a nearly uniform, unique and analytically computable coverage which is in no respect inferior to the one produced by solving the Thomson problem.
To assess the performance of the proposed algorithm under controllable conditions, experiments with simulated data sets have been performed. In this case, the HARDI signals were generated according to model (1) with different values of M (r), D i (r), and s 0 (r) = 1, ∀r. The resulting signals were contaminated by variable levels of Rician noise, giving rise to a set of different SNRs. In this work, we adapt the standard definition of the SNR as SNR = 20 log 10 where s ands denote an original signal and its noise-contaminated version, respectively, and the norms are computed as defined by (18). It is worthwhile noting that the optimal values of regularization parameters λ and µ in (31) are normally a function of the noise level. In the present paper, however, no attempts have been extended to optimize these values for different SNRs. Instead, it was found that λ = 0.03 and µ = 0.05 provided acceptable estimation results in all the simulation scenarios, and hence these values have been used throughout the whole study. Following [37], the scaling parameter ρ in (6) was set to 0.5 and the resolution parameter J in (13) was set to be equal to 1, corresponding to a total of three resolution levels. The number of spherical ridgelet orientations were predefined with m 0 = 4, resulting in M −1 = 16, M 0 = 49 and M 1 = 169 ridgelets spanning the resolution levels j = −1, j = 0 and j = 1, respectively. Thus, the total number of spherical ridgelets used in the reconstruction was equal to 234.
To quantitatively compare the reconstruction results produced by the proposed and references methods for different numbers of sampling directions K and various SNRs, three performance measures were used. The first of the three was the normalized mean-squared error (NMSE) defined as with s(r) being a reference HARDI signal corresponding to location r andŝ(r) being its estimate. Depending on the nature of a specific experiment, the reference signal can be either a simulated signal discretized at 642 spherical points obtained by the 3rd order tessellation of the icosahedron or a signal reconstructed using a maximum possible number of diffusion-encoding orientations. One of the most valuable outcomes of HARDI is in providing an access to computation of orientation distribution functions (ODFs) -the functions whose modes are likely to coincide with the direction of local diffusion flows [14]. Both the SH-based [28] and ridgelet-based [37] methods of reconstruction of HARDI signals come with analytical expressions which relate the HARDI signals to their corresponding ODFs. The latter can in turn be used to recover the directions of local diffusion flows (or, equivalently, the orientations of their related fibre tracts) using, e.g., the steepest ascent procedure detailed in [37]. Suppose u 0 is the true direction of a diffusion flow andũ is its estimate. Then, the angular orientation error δ can be defined (in degrees) as δ = 180 π arccos(u 0 ·ũ).
In this paper, as a performance measure, we use an average angular orientation error which is obtained by averaging the values of δ computed for all "fibres" within a specified Ω d . The last performance measure used in this work is the probability P d of false fibre detection. To define P d , let M (r) be the true number of fibre tracts passing through voxel r (as defined by model (1)). Also, letM (r) be an estimated number of fibres, which is equal to the number of modes (maxima) of the ODF recovered at position r. Then, one can define In addition to the quantitative comparison, the reconstruction results will be evaluated through visual comparison as well. In this paper, we choose to visualize spherical functions by means of 3-D surface plots. Such a plot tends to project away from the origin of R 3 in the directions along which a spherical function is maximized, while passing near the origin in the directions where the function approaches zero.
Finally, our choice of reference methods was motivated by the scope of the main statements made in this paper. First, since we argue that the frame of spherical ridgelets is optimally suited for CS-based reconstruction of HARDI signals, its performance has to be compared with that of alternative representation systems. In particular, the basis of spherical harmonics up to the order 8 inclusive has been used for a different definition of the sensing matrix A in (14). (Note that, in the case of a real and symmetric analysis, this SH-basis consists of 45 functions.) Additionally, the representation system proposed in [46,47] has been exploited in the comparative study as well. This system is formed by applying a set of rotations to a Gaussian kernel of the form d(u) = exp{−b (u T D 0 u)}, with b defined as in (1) and D 0 equal to a (diagonal) diffusion tensor having a mean diffusivity of 766 mm 2 /s and a fractional anisotropy of 0.8. Following [47], the number of rotations (and hence the number of Gaussian basis functions) was set to be equal to 253. For the convenience of referencing, the CS-based reconstruction methods using the spherical ridgelets, the 8th-order spherical harmonics, and the rotated Gaussian kernels will be referred below to as the RDG, SH8 and GSS algorithms, respectively.
To assess the significance of the proposed spatial regularization, all the above algorithms have been applied with two different values of µ in (31), viz. µ = 0 and µ = 0.05. Note that, in the first of these cases, the spatial regularity is ignored, which leads to the sparse-only reconstruction discussed in Section 4.1.
In the second case, on the other hand, the spatial regularity is taken into account and the reconstruction is performed by means of the split Bregman algorithm of Section 5.

In silico experiments
To assess the performance of the proposed and reference methods under controllable conditions, two simulated data sets were used. The first set (referred to below as Phantom #1) had a spatial dimension of 12 × 12 pixels, and consisted of two "fibres" crossing each other at the right angle as it is shown in the upper row of subplots of Fig. 1. In addition, each pixel in the set was assigned an extra diffusion flow in the direction perpendicular to the image plane. As a result, the number of diffusion components M (r) in Phantom #1 varied between 1 and 3. Subsequently, model (1) was used to generate corresponding diffusion-encode images {s k } K k=1 for a range of different values of K. Two different values of b, namely b = 1000 s/mm 2 and b = 3000 s/mm 2 were used for data generation. The diffusion tensors D i (r) in (1) were obtained by applying rotations to a tensor of the form D 0 = diag ([α, β, β]), where α and β were equal to 1700 · 10 −6 and 300 · 10 −6 , respectively. Note that the mean diffusivity and fractional anisotropy of D 0 are equal to 766 mm 2 /s and 0.8, respectively. Thus, the same diffusion tensors were used for data synthesis and for the construction of basis functions in the GSS algorithm, thereby allowing the latter to perform under the best possible conditions.  The lower row of subplots in Fig. 1 depict four examples of the diffusionencoding images obtained for Phantom #1 before their contamination by Rician noise. One can see that the images are piecewise constant functions, which appears to be in a good agreement with the bounded-variation model suggested by (31). However, real images may be more complicated than that. Accordingly, to test the robustness of the proposed regularization scheme, a different in silico phantom was designed. Phantom #2 had a spatial dimension of 16 × 16 pixels and it was obtained through supplementing the configuration of Phantom #1 by an additional circular "fibre" as shown in the upper row of subplots in Fig. 2. The lower row of subplots of the figure show a subset of the resulting diffusionencoded images, which can be seen to no longer exhibit a piecewise constant behaviour characteristic for Phantom #1.
The simulated diffusion-encoded images were contaminated by three different levels of Rician noise, giving rise to SNR of 24, 18 and 12 dB. Some typical examples of the resulting images are demonstrated in Fig. 3, where the upper row of subplots depict a noise-free version of one of the diffusion-encoded images pertaining to Phantom #1 along with its noise-contaminated counterparts. The lower row of subplots in Fig. 3 depict an analogous set of examples for Phantom #2. Observing the figure, one can see that the SNR values have been chosen so as to cover a range of possible noise scenarios, which could be characterized as moderate-to-severe contamination.
As it was mentioned earlier, in our in silico study we compared the performance of three different representation bases, i.e. spherical harmonics (SH8), Gaussian kernels (GSS) and spherical ridgelets (RDG). All the resulting algorithms have been further subdivided into two different types, depending on whether or not the spatial regularization was engaged. Thus, in the absence of the spatial regularization (corresponding to µ = 0), the reconstruction has been performed on a voxel-by-voxel basis, as detailed in Section 4.1. For the convenience of referencing, the corresponding algorithms will be referred to below as SH8-CS, GSS-CS, and RGD-CS. In the case of µ > 0, the estimation has been carried out using the split Bregman method of Section 5. The corresponding algorithms will be referred below as SH8-TV, GSS-TV, and RGD-TV.
The upper subplot of Fig. 4 shows the original field of ODFs of Phantom #1 (corresponding to b = 3000 s/mm 2 ), which have been computed based on Tuch's approximation [14] (i.e. by applying the Funk-Radon transform to the diffusion signals). At the same time, the middle row of subplots of Fig. 4 show the ODFs recovered by (from left to right) SH8-CS, GSS-CS and RDG-CS with K = 16 and SNR = 18 dB. One can see that the inability of the SH basis to sparsely represent HARDI signals results in a poor performance of SH8-CS. A better result is obtained with GSS-CS, which uses a basis of rotated Gaussian kernels, and therefore has a potential to represent the HARDI signals in a sparse manner. Unfortunately, the excessive correlation between the Gaussian basis functions adversely affects the ability of this method to withstand the effect of noise. Consequently, the reconstruction obtained using GSS-CV suffers from sizeable errors. The RDG-CS method, on the other hand, provides an estimation result of a much higher quality, albeit some inaccuracies are still noticeable in the   central part of the phantom. The reconstruction accuracy improves dramatically when the spatial regularization is "switched on", as it is demonstrated by the bottom row of subplots in Fig. 4. Specifically, while SH8-TV is still unable to provide a valuable reconstruction, the estimates obtained using GSS-TV and RDG-TV represent correctly the "flow structure" of Phantom #1. Moreover, among the latter two methods, RDG-TV is clearly the best performer, resulting in a close-to-ideal recovery of the original ODFs. The superiority of RDG-TV over the alternative methods is further evident in the results presented by Fig. 5, which depicts the reconstructions obtained for Phantom #2 (with the same values of b, K and SNR as above).
In general, the reconstruction results obtained using SH8-CS and SH8-TV have been observed to be of a lower quality in comparison to the other methods under consideration. For this reason, in what follows, only the GSS and RDG methods are compared. Thus, Fig. 6 contrasts the performances of GSS-CS, GSS-TV, RDG-CS and RDG-TV in terms of the NMSE criterion. One can see that the best performance here is attained by the RDG-TV algorithm, which results in the smallest values of NMSE for both phantoms and for all the tested values of b, SNR and K. It is also interesting to note that the incorporation of spatial regularization allows GSS-TV to outperform RDG-CS, with the effect of the regularization becoming more pronounced at lower SNRs. On the whole, all the NMSE curves demonstrate an expected behaviour, with the error values increasing proportionally with a decrease in SNR, while going down with an increase in the number of diffusion-encoding gradients K. However, as opposed to the others, the NMSE curves obtained with RDG-TV are characterized by a relatively low rate of convergence, which indicates a reduced sensitivity of RDG-TV to the value of K.
The above algorithms have been also compared in terms of the angular error (42). The results of this comparison are summarized in Fig. 7, which again indicates that the most accurate reconstruction is obtained using the RDG-TV method. In general, the angular error tends to grow as SNR decreases and to converge to a minimum as K increases. As opposed to the case of NMSE, however, there is an additional dependency of the angular error on the type of a phantom in use as well as on the b-value. In particular, the errors obtained for Phantom #2 are (on average) greater than those obtained for Phantom #1. This discrepancy is rooted in the fact that Phantom #2 has a more complex "fibre structure" as compared to Phantom #1. In particular, while the "fibers" of Phantom #1 are designed to cross each other at the right angle, the "fibres" of Phantom #2 are allowed to decussate at much smaller angles, which makes them much harder to resolve. Moreover, this effect becomes more noticeable with a decrease in the b-value, which reduces the resolution of q-ball imaging. Finally, we notice that, on average, GSS-TV performs better than RDG-CS (though still worse than RDG-TV), which justifies the value of spatial regularization.
The comparison in terms of the rate of false fibre detection P d (43) was last in the line of our in silico performance tests; its results are shown in Fig. 8. One can see that, in the case of Phantom #1, RDG-TV yields a virtually zero false detection rate for both values of b, whereas the other methods result in considerably higher values of P d (mainly due to the detection of spurious local maxima in the estimated ODFs). The situation is different for Phantom #2, where all the compared methods yield sizeable errors (especially for b = 1000 s/mm 2 ). However, in comparative terms, the most accurate reconstruction is still obtained by means of the proposed RDG-TV algorithm.

In vivo results
As the next validation step, experiments with real HARDI data were carried out. The proposed algorithm was tested on human brain scans acquired on a 3-Tesla GE system using an echo planar imaging (EPI) diffusion-weighted image sequence. A double echo option was used to suppress eddy-current related distortions. To improve the spatial resolution of EPI, an eight channel coil was used to perform parallel imaging by means of the ASSET technique with a speed-up factor of 2. The data were acquired using 51 gradient directions (quasi-uniformly distributed over the northern hemosphere) with b = 1000 s/mm 2 . In addition, eight baseline (b0) scans were acquired, averaged and used for normalization. The following scanning parameters were used: TR = 17000 ms, TE = 78 ms, FOV = 24 cm, 144 × 144 encoding steps, and 1.7 mm slice thickness. All scans had 85 axial slices parallel to the AC-PC line covering the whole brain.
The main question addressed through the in vivo experiments has been whether or not it is possible to supersede the spatial regularization by prefiltering of HARDI signals. To this end, the RDG-CS algorithm was applied first to the HARDI data containing the full set of K = 51 diffusion gradients. (Note that such dense reconstruction is analogous to the one reported in [37], where the latter is shown to outperform the SH-based estimation [28].) The resulting ODFs have been used as a fiducial against which different reconstruction results were compared.
As the next step, three different subsets of 16, 24 and 32 spherical points were composed out of the original set of 51 diffusion gradients. Within each of these subsets, their corresponding points were chosen so as to result in a quasiuniform coverage of the northern hemisphere. Accordingly, the HARDI data were rearranged into three data sets of size 144 × 144 × 85 × 16, 144 × 144 × 85 × 24 and 144 × 144 × 85 × 32 to emulate compressed sensing data acquisition. The above sets were used to assess the performance of different reconstruction methods. Unfortunately, we have not succeeded to find conditions under which the SH8 and GSS algorithms would provide stable reconstruction results (either with or without pre-filtering). For this reason, only the RDG-CS and RDG-TV algorithms are compared below.
The upper row of subplots in Fig. 9 show the generalized anisotropy (GA) [14] image of a coronal cross-section of the brain along with the reference field of ODFs corresponding to the region indicated by the yellow rectangular. Anatomically, this region is expected to contain the fibre bundles of corona radiata as well as those of superior longitudinal and arcuate fasciculi. The middle row of subplots in the same figure depict the ODFs reconstructed by RDG-CS using    One can see that the quality of reconstruction progressively improves as K increases. It is important to note that, before applying the RDG-CS algorithm, the diffusion-encoded images had been pre-processed by a TV filter to reduce the effect of measurement noises on the estimation result. However, this pre-processing appears to be not nearly as effective as the spatial regularization of the RDG-TV algorithm, whose reconstruction results are shown in the bottom row of subplots in Fig. 9. The above conclusion is further supported by an additional example of Fig. 10, which shows the reconstructions pertaining to the indicated area within an axial cross-section of the brain. (The relevant fibre bundles here are those of cingulum and corpus callosum). As in the previous example, one can see that the most accurate reconstruction is attained by means of the proposed RDG-TV method. The superiority of RDG-TV is also confirmed by the quantitative figures of Table 1, which summarizes the NSME obtained by the compared algorithms for different values of K.

Discussion and Conclusions
When considered as a whole, the HARDI signals which pertain to a given volume of interest can be described as multi-valued (or, more generally, measure-valued) functions from a subset of R 3 to the space of square-integrable spherical functions L 2 (S 2 ). Such functions can be thought of as if they had two "modes of variation" -one in the spatial and another in the diffusion domain. Although applying various inverse problems (a particular instance of which is addressed by the theory of CS) along the spatial and diffusion coordinates independently is not new to the community of medical imaging scientists, formulating a CS reconstruction problem in both domains simultaneously has not been proposed before. Accordingly, the present paper introduced the RDG-TV algorithm which exploits the above idea and can be used for reliable reconstruction of HARDI signals from as few as K = 16 diffusion-encoded scans (as compared to 60-80 scans required by existing reconstruction tools). The algorithm exploits the fact that HARDI signals can be sparsely represented by spherical ridgelets in the diffusion domain, while their associated diffusion-encoded images have bounded variation in the spatial domain. Moreover, it has been shown experimentally that either using different representation bases or excluding the spatial regularization would result in considerably less accurate reconstruction results. At the practical level, the reconstruction is implemented based on the split Bregman approach (with 20 being the maximum number of Bregman iterations used in this study). The resulting algorithm alternates between two estimation stages (38): first, a sequence of basis pursuit de-noising problems are solved independently on a voxel-by-voxel basis, followed by applying a TV filter to a total of K discrete images. Such computations are straightforward to accelerate using a multicore processing, which is another advantage of the proposed reconstruction method. We believe that the algorithm presented in this paper can be improved in a number of ways. First, the square metric used to assess the model fidelity could be replaced by a different metric, which would be more specific to the nature of Rician noise. Second, the fact that diffusion signals are positive-valued could be explicitly incorporated into the reconstruction process in the form of additional constraints. Lastly, the bounded variation model could be substituted by an alternative model, which could (possibly) provide a better account for the spatial regularity of HARDI signals. Exploring the above options constitutes essential part of our ongoing research.
Finally, as the experimental study reported in this paper was comparative in its nature, it was not really important what method to use for approximation of ODFs. Specifically, the present results have been obtained using Tuch's approximation [14]. However, more accurate computation of ODFs is possible based on the solid angle formulation as detailed in [89,90] 3 . It should also be noted that the technique proposed in [89] can be applied to multi-shell HARDI data (i.e. the data acquired for a range of b-values). Until recently, collecting such data has been deemed impractical due to extremely long acquisition time required. We believe, however, that the proposed method for CS-based reconstruction of HARDI data has a potential to help multi-shell HARDI develop into a clinically relevant tool of diagnostic imaging.


Proposal for MWE Annotation in Running Text

We present a proposal for the annotation of multi-word expressions in a 1M corpus of contemporary portuguese. Our aim is to create a resource that allows us to study multi-word expressions (MWEs) in their context. The corpus will be a valuable additional resource next to the already existing MWE lexicon that was based on a much larger corpus of 50M words. In this paper we discuss the problematic cases for annotation and proposed solutions, focus-ing on the variational properties of MWEs.

Introduction
Given the widespread studies of co-occurring words phenomenon, the term 'multi-word expression' (MWE) usually refers to a sequence of words that act as a single unit, embracing all different types of word combinations. Their study is of extreme importance for computational linguistics, where applications find notorious difficulties when dealing with them (Sag et al., 2002).
Having a well-balanced corpus annotated with multi-word expressions offers the possibility to analyze the behavior of MWEs as they appear in running text. Such corpus will contain a rich and diversified set of MWE and also be an excellent resource to evaluate automatic MWE identification systems. Here we propose our approach to the manual annotation of the CINTIL corpus (Barreto et al., 2006) with MWE information. This Portuguese corpus of 1M tokens is a balanced corpus of both spoken and written data from different sources and has been previously annotated with linguistic information such as part-of-speech and lemma and inflection.
As the starting point for our annotation project, we want to use a Portuguese MWE lexicon containing approximately 14,000 entries. The lexicon contains besides idiomatic expressions, also many collocations: expressions of frequently cooccurring words that do not show syntactic or semantic fixedness. We are mostly interested in the idiomatic expressions and will only mark up these in the corpus.

Related Work
There is already quite some work about the creation and representation of MWE lexicons (Baldwin and Kim, 2010). Most of the currently available corpora annotated with MWE information consist of a collection of extracted sentences containing a MWE (for example the data sets in the MWE 2008 shared task 1 ). Fellbaum et al. (2006) report on a larger German example corpus consisting of MWEs with their surrounding sentences. There are also data sets specifically designed for automatic MWE identification, in which part of the sentences contains an idiomatic expression and the other part expresses a literal meaning (e.g. (Sporleder and Li, 2009)). An example of a balanced corpus fully annotated with MWEs is the Prague Treebank which is enriched with a diverse set of MWE annotations (Böhmová et al., 2005).

MWE Lexicon
Our annotation proposal uses information from a lexicon of MWE for Portuguese (available online 2 ). This lexicon is implemented on a MySQL relational database. The MWEs were extracted from a 50M words balanced corpus of Portuguese. The MWE are organized under canonical forms. Also inflectional variations of the canonical forms are recorded, in total the lexicon contains 14,153 canonical forms and 48,154 MWEs variations. For each of those several examples are collected from the corpus. Each MWE entry is also assigned to one or multiple word lemmas, of a total number of 1180 single word lemmas. The MWE were selected from a sorted list of n-grams based on the mutual information measure (Church and Hanks, 1990) and validated manually .

Proposed annotation
In this section we discuss our approach to the annotation of MWEs in the corpus.

Typology
We want to classify each idiomatic MWE occurring in the CINTIL corpus according to a typology that expresses the typical properties of the MWE. Although the lexicon of MWEs covers a wide range of units, from idiomatic expressions to collocations, we decided to restrict our annotation of the corpus to cases of idiomatic MWEs because those are the problematic ones for any task of semantic annotation and disambiguation. The MWE lexicon does not provide labels for idiomatic vs. compositional expressions, so this information will have to be added during the annotation task. Identifying idiomatic MWEs is not a simple task. For clear cases of idiomatic units, the global meaning can not be recovered by the sum of the individual meanings of the elements that compose the expression.
In other cases, only part of the MWE has an idiomatic meaning, while one or more of the elements are used in their literal meaning (e.g saúde de ferro 'iron health'). Deciding if one of the elements of the MWE is literal or not depends in fact of our definition of literal: if we consider it to be the first prototypical meaning of a word, this very restrictive definition will trigger us to label a large number of MWEs as idiomatic. Other MWEs are compositional but receive an additional meaning, like cartão vermelho in football, which is literally a red card but has an additional meaning of punishment.
We want to cover these different cases in our annotation, and to establish a typology that takes into account morpho-syntactic and semantic aspects of the MWE: its functional part-of-speech (PoS) category, the PoS categories of its internal elements, its fixed or semi-fixed nature, its global or partial idiomatic property and motivation, and possible additional meanings.

Division by syntactic category
When studying the MWE lexicon, we noticed different properties of MWEs according to their syntactic patterns. Consequently, we propose to divide our annotation guidelines according to each syntactic pattern and to establish different properties that enables us to distinguish literal from idiomatic usage. At the sentence level, MWEs such as proverbs or aphorisms (e.g.água mole em pedra dura tanto bate até que fura lit. 'water in hard rock beats so long that it finally breaks') have specific properties: they do not accept any possible syntactic changes like passivization or relativization, they do not accept any inflectional variation, the only possible change is lexical (when speakers substitute one or more elements, like we will discuss in section 4.4). However fixed, the meaning of this example is clearly motivated and compositional in the sense that it is recovered by the meaning of the individual elements. On the contrary, MWEs which are verb phrases will admit much more morpho-syntactic variation. Moreover, noun phrases raise specific issues: the most syntactically fixed units will be very close or identical to compound nouns. For example, the meaning of the prepositional modifier of the noun can be literal but the overall expression will still be used as a compound and will denote a very specific entity, frequently from domain-specific languages (projecto de lei 'project of legislation', contrato de compra e venda 'sell contract'). Moreover, the prepositional and adjectival modifiers of the noun will express many different semantic relationships (part of, made of, used for) which interact with the meaning (literal or idiomatic) of the noun (Calzolari et al., 2002). Establishing specific guidelines for these different types of MWEs will enable a more accurate annotation. To decide upon the difficult cases of idiomatic and non-idiomatic usage, we plan to use the intuitions of different annotators.

Linking to MWE lexicon
We will annotate each encountered MWE in the corpus with a link to the MWE-entry in the lexicon, instead of labelling each MWE with its typology. This way we link each MWE to its canonical form and other additional information. Moreover, we can easily gather all occurrences of one particular canonical MWE and check its variation in the corpus. It will also allow us to work with a more detailed typology and will give us the possibility to revise it during the annotation process. It might be difficult to establish beforehand very precise guidelines that will apply to all the MWEs and even to all the MWEs of a specific subtype. Often, guidelines are constantly in need of revision as we encounter slightly different contexts who challenges decisions previously taken. The corpus annotation will enable us to extend the information in the MWE lexicon with typology labels regarding the whole expression (function, idiomatic meaning) but also regarding individual words of the expression as to whether they are obligatory or not.
We plan to add a meaning to idiomatic expressions using a dictionary. We expect that MWEs will be unambiguous: they have the same meaning each time they are used. In some cases, the synonym or paraphrase proposed for the MWE might not be able to replace the MWE in the corpus context. For example, the MWEàs mãos cheias means em grande quantidade 'in large quantity', but this meaning can not always replace the MWE in context.
The annotation process of fully fixed expressions could be retrieved automatically. For the variable expressions we will combine automatic retrieval with manual validation, Here the automatic retrieval step will aim for a high recall and select all sentences that contain the lemmas of the MWE. Without doubt our corpus will contain many MWEs that are not yet listed in the MWE lexicon. Therefore each sentence will need to be checked manually for MWEs. We can create the links between the lexicon and MWEs in the corpus automatically, but again, as not all MWEs will occur in the lexicon, we will need to do a manual validation of the automatic labelling and also add newly discovered MWEs to the lexicon.

MWE Variation
Corpus analysis clearly shows that MWEs have different types of internal variation. Following Moon (1998), we will also assume that, in most of the cases, these expressions "have fixed or canonical forms and that variations are to some extent derivative or deviant". The canonical forms of (variable) expressions are listed in the MWE lexicon. Mapping MWE occurrences in the corpus to their canonical form can be a hard task depending on the flexibility of the MWE. In the next part we discuss our proposal how to handle the annotation of several types of variation in MWEs: lexical, syntactic and structural variation, lexical insertions and truncation of MWEs.

Lexical diversity
MWEs have a wide range of lexical variation and it can apply to any type of grammatical category, although we do notice that verb variation is the commonest type. Studying the lexicon showed us that there is a group of cases in which a word in a MWE can only be replaced by another word from a very limited set (usually not larger than 10 words) of synonyms or antonyms. For these cases this set is already recorded in the MWE lexicon. We mark these variable words as: 'obligatory parts of the MWE and member of a specified list'. In 1 we show an example: the canonical form followed by a sentence containing this MWE and the English translations.
Many MWEs also contain parts that are almost lexically free or only restricted to a semantic class such as person or named entity. These elements are represented in the MWE lexicon with a pronoun (e.g. alguém, algum ('someone', 'something')) or the tag NOUN (with possible gender/number restrictions) when a pronoun cannot substitute the free part. When marking up these elements in the corpus, we will label them with a reference to the pronoun used in the canonical form (example 2).
(1) dizer/ sair da boca para fora (to say / to get out from the mouth outside) Arrependeu-se com o que lhe saiu da boca para fora 'She regretted her slip of the tongue' (2) estar nas mãos de ALGUÉM A nossa vida está nas mãos de Deus 'Our life is in the hands of God' MWEs are not always contiguous: it is frequent to encounter insertion of lexical elements which do not belong to the canonical form of the MWE. Often, the function of the inserted elements is adverbial, quantificational or emphatic. Or the MWE occurs in a negative context, by the insertion of the adverb não. Such inserted elements that are not part of the MWE are not labelled. This is the case of the quantifier muitas in (3), which is not part of the canonical form of the MWE dar voltasà cabeça 'to think'.
(3) Dei muitas voltasà cabeça para encontrar uma solução. 'I've been thinking a lot to find a solution.' Another type of MWE variation is truncation: only a part of the full expression is lexically realized. This phenomenon usually occurs with proverbs and sayings. For example in 4 the bracketed part was not realized in the sentence, but it is part of the canonical form in the MWE lexicon. When marking up such truncated expressions we do not label explicitly this phenomenon, we just mark up the occurring part with a reference link to MWEs in the lexicon.
(4) mais vale um pássaro na mão (do que dois a voar) 'bird in the hand is worth (two in the bush)'

Syntactic variation
An obvious form of syntactic variation is inflection of verbs and nouns. Since Portuguese is a highly inflectional language, practically all the verbs that occur in MWEs inflect, except for some fixed sayings. Also shifting from active to passive voice leads to syntactic variation. We do not label auxiliary verbs as part of the MWE. Several MWEs that have a free part such as example 2 do not only exhibit lexical variation but also syntactic variation: pronominalization (estar nas mãos dele) or with a possessive form (estar nas suas mãos). In such cases we will mark up possessives as part of the MWE but give them an additional label to signal that they are optional elements. However, possessives are not always optional, sometimes it is an obligatory part of the canonical form and we will annotate it normally (e.g. o leão mostra a sua raça. 'the lion shows what he's made off').
Also permutations of the MWE can occur (ex.5). We do not signal this phenomenon in our annotation as this can easily be detected when comparing to the canonical form.

Structural variation
True idioms are both semantically and syntactically fixed. However, language use is creative and can lead to MWEs that only partly match the 'real' MWE as listed in the MWE lexicon. For these cases we mark up the different part with an extra label to clarify which part exactly varies. For example 6.
(6) no pouparé que está o ganho in the saving is the profit no esperar / provar / compararé que está o ganho in waiting / proving / comparing is the profit (7) já dei voltas e voltasà cabeça 'thoughs went on and on in my mind' (8) ALGOé a mãe de todas NOUN-PL 'something is the mother of all x' a educaçãoé a mãe de todas as civilizações a liberdadeé a mãe de todas as virtudes 'education is the mother of all civilizations' 'freedom is the mother of all virtues' Another interesting case is shown in example 7 in wich a part of the MWE is duplicated for emphasis. This should be treated differently than the example in 3. In these cases we will label the duplicated part as 'part of the MWE but optional' (similar to possessives). There are cases in which part of the MWE may vary without any apparent limits, while the other part remains fixed. An example can be found in 8. These are actually just an extension of ones we already discussed (see example 2) and we treat them in the same matter.

Conclusion
In sum, we propose to split the annotation of MWEs to develop separate annotation guidelines for the grammatical categories, as we have observed that e.g. nominal MWEs behave differently than verbal MWEs. Each MWE in the running text will be linked to its canonical form in the lexicon. The lexicon itself will be enhanced with additional information such as typology information and MWE meaning. Special elements of the MWE such as optional or variable parts will be explicitly marked as such both in the lexicon and in the annotation of the MWE in the corpus. We are convinced that the implementation of our proposal will lead to a rich new resource that can help us study the behavior of MWE in more depth. We also plan to use this resource for the development and evaluation of automatic MWE identification systems.