Eugen Barbu
email: eugen.barbu@univ-rouen.fr

Pierre Héroux

Sébastien Adam

Éric Trupin

Clustering document images using a bag of symbols representation

Document image classification is an important step in document image analysis. Based on classification results we can tackle other tasks such as indexation, understanding or navigation in document collections. Using a document representation and an unsupervised classification method, we may group documents that from the user point of view constitute valid clusters. The semantic gap between a domain independent document representation and the user implicit representation can lead to unsatisfactory results.

In this paper we describe document images based on frequent occurring symbols. This document description is created in an unsupervised manner and can be related to the domain knowledge. Using data mining techniques applied to a graph based document representation we find frequent and maximal subgraphs. For each document image, we construct a bag containing the frequent subgraphs found in it. This bag of "symbols" represents the description of a document.

We present results obtained on a corpus of 60 graphical document images.

Introduction

A document image analysis (DIA) system transforms a document image into a description of the set of objects that constitutes the information on the document and which are in a format that can be further processed and interpreted by a computer [START_REF] Antonacopoulos | Introduction to Document Image Analysis[END_REF]. Documents can be classified in mostly graphical or mostly textual documents [START_REF] Nagy | Twenty years of document analysis in pami[END_REF]. The mostly textual documents also known as structured documents respect a certain layout and powerful relations exist between components. Examples of such documents are technical papers, simple text, newspapers, program, listing, forms,. . . Mostly graphical documents do not have strong layout restrictions but usually relations exist between different document parts.

Examples of this type of documents are maps, electronic schemas, architectural plans. . . For these two categories of documents, graph based representations can be used to describe the image content (e.g. region adjacency graph [START_REF] Pavlidis | Algorithms for Graphics and Image Processing[END_REF] for graphical and Voronoibased neighborhood graph [START_REF] Bagdanov | Fine-grained document genre classification using first order random graphs[END_REF] for textual document images).

In this paper we present an approach similar with the "bag of words" method from Information Retrieval (IR) field. We describe a document using a bag of symbols found automatically using graph mining [START_REF] Washio | State of the art of graph-based data mining[END_REF] techniques. In other words, we consider the frequent subgraphs of a graph-based document representation as "symbols" and we investigate whether the description of a document as a bag of "symbols" can be profitably used in a clustering task.

The approach has the ability to process document images without knowledge of, or models for, document content. In the literature one can find papers dealing with clustering of textual documents using frequent items [START_REF] Fung | Hierarchical document clustering using frequent items[END_REF] and description of XML documents using frequent trees [START_REF] Termier | Mining xml data with frequent trees[END_REF] but we do not know any similar approach in the DIA field.

The motivation for our study is the fact that unsupervised classification can represent the starting point for semisupervised classification or indexation and retrieval from document collections. Also, the existing clustering solutions for document images are usually domain dependent and can not be used in an "incoming document flow" (fax, business mail. . . ) setting, where supervised techniques are not at hand. The outline of this paper is as follows. In section 2 we present a graph representation and how we create this representation from a document image. Section 3 presents the graph-mining method used, in section 4 we describe how we create clusters based on dissimilarities between bags of symbols. Section 5 shows some experimental results. We conclude the paper and outline perspectives in section 6.

Document graph based representations

Eight levels of representation in document images are proposed in [START_REF] Blostein | Document representations[END_REF]. These levels are ordered in accordance with their aggregation relations. Data array level, primitive, lexical, primitive region, functional region, page, document, and corpus level are the representation levels proposed.

Without loosing generality, in the following paragraphs we focus our attention on a graph-based representation build from the primitive level. The primitive level contains objects such as connected components (set of adjacent pixels with the same color) and the relations between them.

Let I be an image and C(I) the connected components from I, c ∈ C(I) is described as c = (id, P ), where id is a unique identifier and P the set of pixels the component contains. Based on this set P , we can compute the center for the connected component bounding box and also we can associate a feature vector to it. Considering that, c = (id, x, y, v), v ∈ R n . Subsequently using a clustering procedure on the feature vectors we can label the connected component and reach the description C = (id, x, y, l) where l is a nominal label. The graph G(I) representing the image is G = G(V (I), E(I)). Vertices V (I) correspond to connected components and are labeled with component labels. An edge (u, w) between vertex u and vertex w exists if f (u.xw.x) 2 + (u.yw.y) 2 < t, where t is a threshold that depends on the image I global characteristics (size, number of connected components. . . ).

The exact methodology employed to construct the graph representation is subsequently presented. From a binary (black and white) document image we extract connected components (black and white). The connected components will be the graph nodes. For each connected component we extract features. In the actual implementation the extracted characteristics are rotation and translation invariant features based on Zernike moments [START_REF] Khotazad | Invariant image recognition by zernike moments[END_REF]. The invariants represent the magnitudes of a set of orthogonal complex moments of a normalized image.

The following step is to associate each connected component a label.

The two main categories of clustering methods are partitional and hierarchical. Partitional methods can deal with large sets of objects ("small" in this context means less than 300) but needs the expected number of clusters in input. Hierarchical methods can overcome the problem of number of clusters by using a stopping criterion [START_REF] Gordon | Classification. Chapman & Hall[END_REF] but are not applicable on large sets due to their time and memory consumption. In our case the number of connected components that are to be labeled can be larger than the limit of applicability for hierarchical clustering methods . In the same time we cannot use a partitional method because we do not know the expected number of clusters. Based on the hypothesis that a "small" sample can be informative for the geometry of data, we obtain in a first step an estimation for the number of clusters in data. This estimation is made using an ascendant clustering algorithm with a stopping criterion. The number of clusters found in the sample is used as input for a partitional algorithm applied on all data. We tested this "number of cluster estimation" approach using a hierarchical ascendant clustering algorithm [START_REF] Gordon | Classification. Chapman & Hall[END_REF] that uses Euclidean distance to compute the dissimilarity matrix, complete-linkage to compute between-clusters distances, and Calinsky-Harabasz index [START_REF] Milligan | An examination of procedures for determining the number of clusters in a data set[END_REF] as a stopping criterion. The datasets (T 1 , T 2 , T 3 ) 2 we present predicted cluster numbers obtained for different sample sizes. If the test set is T and |S| = 50, after repeating ten times the sampling procedure we obtain a set of estimations for the number of clusters. We can see that by using a majority voting decision rule we can find the good number of clusters in most of the cases and even when the sample size is very small (50 or 100) compared with the data set size.

(
We employed our sampling approach combined with the k-medoids clustering algorithm [START_REF] Kaufmann | Statistical Data Analysis based on the L1 Norm and Related Methods, chapter Clustering by Means of Medoids[END_REF] on the connected components data set from images in our corpus (see section 5). The k-medoids clustering algorithm is a more robust version of the well known kmeans algorithm. The images from our corpus contains 6730 connected components. The proposed number of clusters using ten samples of size 600 is [START_REF] Washio | State of the art of graph-based data mining[END_REF][START_REF] Tabbone | Matching of graphical symbols in line-drawing images using angular signature information[END_REF]17,[START_REF] Washio | State of the art of graph-based data mining[END_REF][START_REF] Washio | State of the art of graph-based data mining[END_REF]19,[START_REF] Gordon | Classification. Chapman & Hall[END_REF]17,[START_REF] Termier | Mining xml data with frequent trees[END_REF][START_REF] Washio | State of the art of graph-based data mining[END_REF] and by considering the majority we use 16 clusters as input to the partitional clustering algorithm.

After labeling the connected components (nodes in the graph) subsequently we describe the way we add edges to the graph. The edges can be labeled or not (if unlabeled the significance is Boolean : we have or not a relation between two connected components) and can be relations of spatial proximity, based on "forces" [START_REF] Tabbone | Matching of graphical symbols in line-drawing images using angular signature information[END_REF], orientation or another criterion. In our actual implementation the distance between centers of connected components is used (see Fig. 1). If the distance between two connected components centers is smaller than a threshold, then an edge will link the two components (nodes).

Table 2. Proposed cluster numbers

Figure 1. An image and its associated graph

Graph Mining

"The main objective of graph mining is to provide new principles and efficient algorithms to mine topological substructures embedded in graph data" [START_REF] Washio | State of the art of graph-based data mining[END_REF].

Mining frequent patterns in graphs is the problem of finding in a set of graphs those subgraphs that occur more times than a threshold (minimum support). Because the number of patterns can be exponential this problem complexity can also be exponential. An approach to solve this problem is to start with finding all frequent patterns with one element, then all patterns with two elements, etc. in a level-by-level setting . In order to reduce the complexity different constraints are used: the minimum support, the subgraphs are connected, and not overlapped.

An important concept is that of maximal subgraph. A graph is said to be maximal if it does not have a frequent super-graph.

In our document image analysis context we are interested in finding maximal frequent subgraphs because we want to find symbols but to ignore their parts.

A system that is used to find frequent patterns in graphs is FSG (Frequent Subgraph Discovery) that "finds patterns corresponding to connected undirected subgraphs in an undirected graph database" [START_REF] Seno | Pafi, a pattern finding toolkit[END_REF].

The input for the FSG program is a list of graphs and a minimum support threshold. Each graph represents a transaction. We present subsequently how we construct the transaction list starting from a set of document images. Using the procedure presented in section 2 we create for each document an undirected labeled graph. Every connected component of this graph represents a transaction. Using FSG we extract the frequent subgraphs and we construct a bag of frequent subgraphs occurring in each document.

Figure 2. Occurences of a frequent subgraph in an image

In the following paragraphs we consider that the frequency condition is sufficient for a group of connected components to form a symbol and we will conventionally make an equivalence between the frequent subgraphs found and symbols. As we can see in the example (Fig. 2) the proposed symbols are far from being perfect due to the image noise, connected components clustering procedure imperfections,. . . however we can remark the correlation between this artificial symbol and the domain symbols.

Dissimilarity between document descriptions

A collection of documents is represented by a symbolby-document matrix A, where each entry represents the occurrences of a symbol in a document image, A = (a ik ) where a ik is the weight of symbol i in document k. Let f ik be the number of occurrences of symbol i in document k, N the number of documents in the collection, and n i the total number of times symbol i occurs in the whole collection. In this setting conform with [START_REF] Dumais | Improving the retrieval information from external ressources, behaviour research methods[END_REF] one of the most effective weighting scheme is entropy-weighting. The weight for symbol i in document k is given by :

a ik = log (1 + f ik ).   1 + 1 log N n j=1 f ij n i log f ij n i  
Now, considering two documents A, B with the associated weights A = (a 1 , a 2 , . . . , a t ), B = (b 1 , b 2 , . . . , b t ) where t is the total number of symbols, then

d(A, B) = 1 - t i=1 a i .b i t i=1 a 2 i . t i=1 b 2 i
represents a dissimilarity measure based on the cosine correlation.

Experiments

A comparison between results obtained using the proposed document representation and three other representations is made in the following paragraphs. On a corpus of graphical document images we have extracted different sets of features. Each document image is described with one of the following types of features : Zernike moments for the whole image (a vector with 16 components), pixel densities (the feature vector considered is composed of the 85 (1+4+16+64) gray levels of a 4-level-resolution pyramid [START_REF] Ballard | Computer Vision[END_REF], see Fig 3), connected components label list, and symbol label list.

(i, j) = |C j ∩ C ′ i | C ′ i the precision prec(i, j) = |C j ∩ C ′ i | C j and F ij =
2.prec(i, j).rec(i, j) prec(i, j) + rec(i, j)

represents the F-Measure. The overall F-Measure of a clustering is

F = l i=1 |C ′ i | |D| max j=1...k F ij
F-measure is 1.0 if the matching between the two partitions (ground truth and the one proposed by the clustering algorithm) is perfect.

Our corpus contains 30 images from the class of a French telephony operator (FT) maps, 25 electronic schemas, and 5 architectural plans. These images are scanned images that contain real and artificial noise.

We can see that the connected component list approach obtains good results compared with the simple approaches

Table 3. Experimental results

(Zernike moments and densities). In the same time the symbol list approach representation is more compact than the connected component list and also obtains better results.

Conclusions

The research undertaken represents a novel approach for clustering document images. The approach uses data mining tools for knowledge extraction. It automatically finds frequent symbols. These frequent patterns are part of the document model and can be put in relation with the domain knowledge. The exposed method can be applied to other graph representations of a document. In the near future, we will apply this approach to layout structures of textual document images. Another follow up activity is to quantify the

Figure 3 .

 3 Figure 3. Four level resolution pyramid

Figure 4 .

 4 Figure 4. Images from the corpus

Table 1 . Data sets description

 1 see Table1) are synthetically generated and contains well separated (not necessary convex) clusters.

	T	|T |	number of clusters
	T 1 24830	5
	T 2 32882	15
	T 3 37346	24

Considering S the sample extracted at random from a test set, in Table

Table 4 .

 4 Confusion matrix detailsway noise can affect the connected components labeling, and the manner in which an incorrect number of clusters can affect the graph mining procedure. Based on this error propagation study we can ameliorate our method. Other possible improvements can be obtained if we would employ a graph-based technique that can deal with error tolerant graph matching.


Michael Sebbagh
email: <michael.sebbagh@inserm.fr

Sylviane Olschwang

Santoni

Jean-Paul Borg

The LKB1 complex-AMPK pathway: the tree that hides the forest

Keywords: Humans, Neoplasms, physiopathology, Protein Kinases, genetics, metabolism, Protein-Serine-Threonine Kinases, genetics, metabolism

Initially identified as the PAR-4 homologue, the serine threonine kinase LKB1 is conserved throughout Caenorhabditis elegans evolution and ubiquitously expressed. In humans, LKB1 is causally linked to the Peutz-Jeghers syndrome (PJS) and is one of the most commonly mutated genes in several cancers like lung and cervical carcinomas. These observations have led to classify as tumour lkb1 suppressor gene. Although considerable dark zones remain, an impressive leap in the understanding of LKB1 functions has been done during the last decade. Role of LKB1 as a major actor of the AMPK/mTOR pathway connecting cellular metabolism, cell growth and tumorigenesis has been extensively studied probably to the detriment of other functions of equal importance. This review will discuss about LKB1 activity regulation, its effectors and clues on their involvement in cell polarity.

Introduction

Initially identified as the PAR-4 homologue involved in early asymmetric cell division , the vertebrate Caenorhabditis elegans [START_REF] Kemphues | Identification of genes required for cytoplasmic localization in early C. elegans embryos[END_REF] LKB1 gene , also called XEEK1 in , is conserved throughout evolution from worms to mammals. In humans, STK11 [START_REF] Jenne | Peutz-Jeghers syndrome is caused by mutations in a novel serine threonine kinase[END_REF] Xenopus [START_REF] Su | Cloning and characterization of a novel serine/threonine protein kinase expressed in early Xenopus embryos[END_REF] gene locates on chromosome 9p13.3 and includes 11 exons which encode for two LKB isoforms by alternative splicing , . STK11 [START_REF] Jenne | Peutz-Jeghers syndrome is caused by mutations in a novel serine threonine kinase[END_REF] [4 5 ] The long LKB1 (50kDa) form is detectable from 7 days old embryo and is ubiquitously expressed in adult with a notable higher [START_REF] Rowan | In situ analysis of LKB1/STK11 mRNA expression in human normal tissues and tumours[END_REF] expression in pancreas, liver and skeletal muscle , . The shorter form (48kDa), which lacks the last 29 C-terminus residues, appears [6 7 ] restricted to testis . Interest in LKB1 has dramatically increased when it was found to be causally linked to the Peutz-Jeghers syndrome [START_REF] Towler | A novel short splice variant of the tumour suppressor LKB1 is required for spermiogenesis[END_REF] , , a rare autosomal dominant disorder characterized by hamartomatous polyposis and mucocutaneous melanin pigmentation .

[2 8 ] [9 ]
Indeed, in 70 of PJS cases , LKB1 kinase activity appears lost due to genomic deletions or single mutations affecting either % [10 ] [START_REF] Launonen | Mutations in the human LKB1/STK11 gene[END_REF] transcript splicing or protein primary structure , . Patients affected by this syndrome have a 18 times higher relative risk of [START_REF] Hastings | An LKB1 AT-AC intron mutation causes Peutz-Jeghers syndrome via splicing at noncanonical cryptic splice sites[END_REF] [ 11 13 ] developing cancers at multiple sites, mostly lung, pancreas, ovary, breast and colon than the general population , . This had led to [14 15 ] define whether sporadic mutations of gene could occur in tumours. At the exception of non-small-cell lung cancer where is LKB1 LKB1

found mutated in 33 of cases, somatic mutations of this gene are rather rare . Nevertheless, its haploinsufficiency has been clearly % [START_REF] Sanchez-Cespedes | A role for LKB1 gene in human cancer beyond the Peutz-Jeghers syndrome[END_REF] involved in poor outcome for pancreas , breast , endometrial and liver adenocarcinomas . These observations have led to [START_REF] Morton | LKB1 haploinsufficiency cooperates with Kras to promote pancreatic cancer through suppression of p21-dependent growth arrest[END_REF] [18 ]

[19 ] [START_REF] Kim | Genetic analysis of the LKB1/STK11 gene in hepatocellular carcinomas[END_REF] classify LKB1 as a tumour suppressor which is a rare feature for a kinase . [START_REF] Sanchez-Cespedes | A role for LKB1 gene in human cancer beyond the Peutz-Jeghers syndrome[END_REF] Generation of LKB1 knock out (KO) mouse models was the most informative to determine LKB1 functions. Thus, homozygous LKB1 KO mice die at 10.5 days of development mainly due to defect in neural tube closure and vascular abnormalities attributed, for this latter, to a VEGF level expression . At the heterozygous state, mice recapitulate PJS with an average tumour free survival of 43 weeks [START_REF] Ylikorkala | Vascular abnormalities and deregulation of VEGF in Lkb1-deficient mice[END_REF] . In addition, conditional LKB1 KO mouse models targeting numerous organs confirmed LKB1 suppressor activity in breast , [22 ] [23 ] pancreas , prostate and skin where tumours appearance is correlated with an altered cell organization such as the nucleus [START_REF] Hezel | Pancreatic lkb1 deletion leads to acinar polarity defects and cystic neoplasms[END_REF] [25 ] [START_REF] Gurumurthy | LKB1 deficiency sensitizes mice to carcinogen-induced tumorigenesis[END_REF] positioning in Langerhans -cells in the case of pancreas , . In the same manner, LKB1 deficiency in the thymus prevents mainly through phosphorylation, at a conserved N-terminal residue, of 14 members of the AMPK family kinase which thus become catalytically actives , so driving this wide range of physiological processes. In recent years, research on LKB1, especially on its regulation mode as well as on its driven functions at cellular level, has been greatly intensified in order to define how this kinase exerts its tumour suppressor properties. This review focuses on current knowledge about LKB1 activity regulation, its effectors and clues on their involvement in cell polarity.

Active LKB1 is a complex

Initial investigations on LKB1 functions and regulation mechanisms were difficult because of the poor LKB1 intrinsic kinase activity.

However, identification of the pseudo kinase STRAD (STe-20 Related ADaptor) as a LKB1 partner was a real breakthrough since this association dramatically improves LKB1 kinase activity . Interaction between STRAD and LKB1 appears to be stabilized by a third [START_REF] Baas | Activation of the tumour suppressor kinase LKB1 by the STE20-like pseudokinase STRAD[END_REF] protein, MO25 (mouse protein25) . This has led to establish that, in fact, active LKB1 should be considered as a heterotrimeric [START_REF] Boudeau | MO25alpha/beta interact with STRADalpha/beta enhancing their ability to bind, activate and localize LKB1 in the cytoplasm[END_REF] complex (1:1:1) hereafter referred to as LKB1 complex ( ). figure 1 In lower organisms where the STRAD and MO25 homologues are expressed, even though this has not been formally established yet, the LKB1 complex formation should be conserved. Indeed in those organisms, MO25 and STRAD homologues contain the critical residues required for their activity such as the STRAD C-terminal WEF motif . In mammals, an additional level of complexity occurs.

[44 ]

Indeed, two human STRAD and MO25 paralogues, defined for both as and , have been characterised and appear ubiquitously expressed α β

. Furthermore, at least 4 STRAD and 2 STRAD isoforms all derived by alternative splicing have been defined. Therefore, [START_REF] Boudeau | MO25alpha/beta interact with STRADalpha/beta enhancing their ability to bind, activate and localize LKB1 in the cytoplasm[END_REF] α [START_REF] Marignani | Novel splice isoforms of STRADalpha differentially affect LKB1 activity, complex assembly and subcellular localization[END_REF] β [START_REF] Sanna | ILPIP, a novel anti-apoptotic protein that enhances XIAP-mediated activation of JNK1 and protection against apoptosis[END_REF] cells may express more than one kind of LKB1 complex and in theory could reach up to 16 different ones. Interestingly, some C-terminus STRAD splice variants, missing their MO25 binding motif, keep their ability to induce LKB1 kinase activity . Nevertheless, α in vitro [START_REF] Marignani | Novel splice isoforms of STRADalpha differentially affect LKB1 activity, complex assembly and subcellular localization[END_REF] this aspect remains unexplored, as LKB1 studies focus almost exclusively on the LKB1-STRAD -MO25 complex, but should be kept in

α α
mind since all these potential complexes could have specific functions or regulation modes.

The structure of the LKB1-STRAD -MO25 complex has been recently resolved and reveals that LKB1 kinase activity results from

α α
an allosteric modification induced by STRAD . Although unable to have an intrinsic kinase activity, due to absence of critical α [START_REF] Zeqiraj | Structure of the LKB1-STRAD-MO25 complex reveals an allosteric mechanism of kinase activation[END_REF] (Asp-Gly-Phe) residues to coordinate magnesium binding in its subdomain VII , STRAD adopts a conformation close to an active [START_REF] Baas | Activation of the tumour suppressor kinase LKB1 by the STE20-like pseudokinase STRAD[END_REF] α " "

kinase which associates with LKB1 as one of its substrates . In addition, this reveals an unrecognized STRAD role by promoting [START_REF] Zeqiraj | Structure of the LKB1-STRAD-MO25 complex reveals an allosteric mechanism of kinase activation[END_REF] α interaction between MO25 and LKB1. Indeed, MO25 plays a crucial role in stabilizing the conformation of the LKB1 activation loop α α required for the kinase activity induced by STRAD binding . Thus, this explains how lack of the LKB1 phosphorylation site in the α [START_REF] Zeqiraj | Structure of the LKB1-STRAD-MO25 complex reveals an allosteric mechanism of kinase activation[END_REF] activation loop, generally required to induce and stabilize the switch from inactive to active kinase, is dispensable. This feature leads to the conclusion that LKB1 needs to be associated with STRAD and MO25 to be active since no LKB1 autophosphorylation would be able to α α maintain this active conformation .

[44 ]

Regulation of LKB1 complex

Although few things have been clearly established on this point, the fact that LKB1 can phosphorylate and activate almost all known members of the AMPK family kinase implies that the LKB1 complex should be sharply regulated to trigger specific responses. Crystal structure had confirmed that LKB1 complex activity is regulated through the LKB1/STRAD interaction. Since expression level or stability of STRAD or LKB1 proteins does not appear to be involved in LKB1 complex regulation, the eventuality that post-translational modifications could be involved has been investigated.

Post-translational modifications

Several phosphorylation sites have been defined on LKB1. Four threonines phosphorylated through LKB1-STRAD interaction are considered as autophosphorylation sites revealing LKB1 complex activity . In the same conditions, STRAD phosphorylation sites [START_REF] Baas | Activation of the tumour suppressor kinase LKB1 by the STE20-like pseudokinase STRAD[END_REF] α have also been identified, though these sites are not conserved in the STRAD paralogue . LKB1 Thr363 is phosphorylated after DNA β [START_REF] Zeqiraj | Structure of the LKB1-STRAD-MO25 complex reveals an allosteric mechanism of kinase activation[END_REF] damage in an ATM kinase-dependent manner , . This phosphorylation does not modify LKB1 complex catalytic activity, but has [47 48 ] been recently demonstrated to be required for inactivation of CRTC2, a transcriptional coactivator of CREB, a critical component of germinal center for B cell proliferation . Even if LKB1 Ser325 and 428, respectively phosphorylated by at least ERK and p90RSK, do [START_REF] Sherman | AID-induced genotoxic stress promotes B cell differentiation in the germinal center via ATM and LKB1 signaling[END_REF] not affect LKB1 complex kinase activity , they have interestingly been described to impede LKB1 inhibitory effect in [49 50 ] anchorage-independent growth, concomitantly to a reduction of AMPK activation . Phosphorylation on several other LKB1 sites [49 -51 ] like Ser31 or Thr189 has also been observed, but kinases involved and functional consequences are still unknown . However, these [START_REF] Sapkota | Identification and characterization of four novel phosphorylation sites (Ser31, Ser325, Thr336 and Thr366) on LKB1/STK11, the protein kinase mutated in Peutz-Jeghers cancer syndrome[END_REF] sites do not appear to be involved in LKB1 complex activity since their substitution by alanine does not significantly affect the LKB1 kinase activity measured . LKB1 can also be farnesylated at a C-terminal CAAX sequence only present in the LKB1 long in vitro [START_REF] Sapkota | Identification and characterization of four novel phosphorylation sites (Ser31, Ser325, Thr336 and Thr366) on LKB1/STK11, the protein kinase mutated in Peutz-Jeghers cancer syndrome[END_REF] form. This modification gives rise to LKB1 higher affinity for membrane localization, as alanine replacement of the critical cysteine in CAAX sequence abolishes it , , . Nevertheless, farnesylation does not affect LKB1 complex activity , . The only LKB1 [47 53 54 ] [49 51 ] post translational modification referred to modulate kinase activity of the complex is acetylation. Indeed, it has been reported that LKB1 could be acetylated on nine different lysines. Among those, acetylation of the LKB1 Lys48 reduces its affinity for STRAD and then α decreases activity of the complex . Although the acetyl transferase responsible for this modification has not been identified yet, SIRT1 [START_REF] Lan | SIRT1 modulation of the acetylation status, cytosolic localization, and activity of LKB1. Possible role in AMP-activated protein kinase activation[END_REF] appears to be as the deacetylase able to remove this acetyl group and then favour LKB1 complex formation and activity . These [START_REF] Lan | SIRT1 modulation of the acetylation status, cytosolic localization, and activity of LKB1. Possible role in AMP-activated protein kinase activation[END_REF] observations suggest that LKB1 complex activity is not mainly regulated by a post translational modification, even though some of them have functional consequences . Furthermore, to our knowledge, no studies have still observed LKB1 complex activation under any [START_REF] Zheng | Oncogenic B-RAF negatively regulates the tumor suppressor LKB1 to promote melanoma cell proliferation[END_REF] stimulating cues, leading to the hypothesis that LKB1 complex could be constitutively active , and that its regulation may [51 53 ] dependent on its intracellular localization .

[53 ]

Subcellular localisations

LKB1 complex localization in cells was initially investigated in over-expressing systems using tagged constructs , . It was [42 56 -58 ] then observed that LKB1 predominantly localizes in the nucleus and to lesser extent in the cytosol. A nuclear localization signal (NLS) at the N-amino terminal of LKB1 was found to be critical for LKB1 nuclearization . Interestingly, overexpressed STRAD is clearly [56 -58 ] α able to displace LKB1 from the nucleus to the cytoplasm compartment , , whereas STRAD was unable to do so . Thus, it [42 43 56 ] β [START_REF] Dorfman | STRADalpha regulates LKB1 localization by blocking access to importin-alpha, and by association with Crm1 and exportin-7[END_REF] was proposed that active LKB1/STRAD was cytosolic whereas the nuclear LKB1 was active through its interaction with STRAD .

α β [START_REF] Dorfman | STRADalpha regulates LKB1 localization by blocking access to importin-alpha, and by association with Crm1 and exportin-7[END_REF] Nevertheless, recent results challenge this model. Indeed, improvement of LKB1 serological tools allowed to immunostain endogenous LKB1. Results from those experiments, confirmed by cell fractionation assays, show that endogenous LKB1 localizes mainly in the cytosolic and membrane fractions and appears to be absent from the nucleus , . Failure to detect endogenous nuclear LKB1 has been [4 53 ] reported in different cell types, leading to reconsider the functional relevance of the identified NLS sequence. In addition, while LKB1 functions appear to be largely conserved throughout evolution, its or homologues do not have NLS sequence and C. elegans Drosophila are not significantly found in the nucleus .

[54 ]

Endogenous LKB1 immunostaining in migrating cells shows LKB1 concentration at the leading edge . In the same way, in fully [START_REF] Zhang | The tumor suppressor LKB1 regulates lung cancer cell polarity by mediating cdc42 recruitment and activity[END_REF] polarized plasma membrane has also been observed in . In , although STRAD has not been localized yet, LKB1 was C.elegans [START_REF] Narbonne | Differential requirements for STRAD in LKB1-dependent functions in C. elegans[END_REF] Drosophila observed at the cell cortex , . Functionally, the LKB1 complex presence at the basolateral domain is correlated to its ability to [54 60 ] activate one of its substrate, AMPK , suggesting that LKB1 complex regulation could be governed through its intracellular localization [START_REF] Sebbagh | Regulation of LKB1/STRAD localization and function by E-cadherin[END_REF] allowing proximity with its substrates. This eventuality is strengthened by sparse observations in which LKB1 Ser431 phosphorylation reduces its affinity for the membrane as well as its ability to properly activates AMPK .

[4 ] [START_REF] Zheng | Oncogenic B-RAF negatively regulates the tumor suppressor LKB1 to promote melanoma cell proliferation[END_REF] Altogether, this gives rise to the idea that the LKB1 complex is constitutively active and regulation of cellular processes in which it is involved appears dependent of its subcellular localizations. Very recently, a second LKB1 localization identified in primary cilium of polarized epithelial cells strengthen this hypothesis, nevertheless mechanism(s) regulating these spatial segregations remain to be [START_REF] Boehlke | Primary cilia regulate mTORC1 activity and cell size through Lkb1[END_REF] defined.

LKB1 complex effectors and substrates

Effectors

LKB1 effectors have been identified mainly by mass spectrometry and yeast two hybrid assays. Thus, LKB1 was found to interact with the Hsp90/Cdc37 protein complex contributing likely to its cellular stability , . LKB1 was also reported to interact with LIP1 [62 63 ] (LKB1 interacting protein 1) which contains 6 Leucine rich repeats (LRR) . Function of LIP1 and consequences of its interaction with [START_REF] Smith | LIP1, a cytoplasmic protein functionally linked to the Peutz-Jeghers syndrome kinase LKB1[END_REF] LKB1 remain elusive even though effects on TGF signalling pathway have been suggested . In the same manner, LKB1 interaction β [START_REF] Smith | LIP1, a cytoplasmic protein functionally linked to the Peutz-Jeghers syndrome kinase LKB1[END_REF] with PTEN was observed without clear evidence of functional consequences. Although poorly investigated in the LKB1 context, [START_REF] Mehenni | LKB1 interacts with and phosphorylates PTEN: a functional link between two proteins involved in cancer predisposing syndromes[END_REF] potentially through LIP1, TGF and PTEN have been described to be involved in several syndromes like juvenile polyposis which leads to β formation of gastrointestinal hamartomas bearing different histological features . Nevertheless, since mechanisms responsible for the [START_REF] Merg | Genetic conditions associated with intestinal juvenile polyps[END_REF] formation of those benign polyps are still undefined, deregulation of the cross talk between TGF , through SMAD4, and mTOR pathway β negatively regulated by PTEN and LKB1 , could play a role in these pathologies. In addition, LKB1 was also described to interact [START_REF] Shackelford | The LKB1-AMPK pathway: metabolism and growth control in tumour suppression[END_REF] with Estrogen receptor-(ER ) and with Brahma related gene-1 (BRG1) , both involved in transcriptional regulation.

α α [67 ] [68 ]
Nevertheless, LKB1 does not appear to phosphorylate these latest partners suggesting perhaps potential indirect LKB1 effect , . [67 68 ] LKB1 complex substrates, the AMPK related kinase family (AMP kinase) was the first LKB1 substrate identified and it is, by far, the one to have been better characterised , . This AMPK [69 70 ] kinase is a heterotrimer composed of a catalytic (AMPK ) and two regulatory (AMPK and AMPK ) subunits. AMPK is activated when

α β γ
the intracellular AMP/ATP ratio increases, for instance during hypoxia or nutrient deprivation. In these conditions, LKB1 complex kinase activity does not increase but cytosolic accumulated AMP binds to AMPK subunit allowing, by allosteric modification, the γ [START_REF] Scott | CBS domains form energy-sensing modules whose binding of adenosine ligands is disrupted by disease mutations[END_REF] accessibility to the critical threonine of the AMPK T-loop kinase domain targeted by LKB1 complex, so leading to its phosphorylation and then to AMPK activation. Main interest in this LKB1 substrate comes from the indication that AMPK modulates the activity of multiple downstream targets such as Raptor (Regulatory associated protein of mTOR) which, in turn, repress the mTORC1 pathway . mTORC1 [START_REF] Shackelford | The LKB1-AMPK pathway: metabolism and growth control in tumour suppression[END_REF] integrates growth factors and nutriment inputs and controls eukaryote cell growth deregulated in most human tumours . Thus, LKB1 [START_REF] Guertin | Defining the role of mTOR in cancer[END_REF] could be the link between energy metabolism and tumour growth suspected since more than 50 years . This hypothesis is particularly [START_REF] Fogarty | C-terminal phosphorylation of LKB1 is not required for regulation of AMP-activated protein kinase, BRSK1, BRSK2, or cell cycle arrest[END_REF] attractive to explain the mechanism by which LKB1 exerts its tumours suppressor activity. However, AMPK activation could occur by multiples alternative ways than through the LKB1 complex. Indeed, AMPK can be activated by alternative kinases such as CAMKKβ (Ca2 /calmodulin-dependent protein kinase kinase ) or TGF -activated kinase-1 (TAK1) which probably contribute to some

+ β [73 ] β [74 ]
consequences attributed to LKB1 complex. Details about AMPK regulated pathways, especially its involvement in mTOR pathway and cell contractility regulation, are reviewed by Tomi Makela and al. in this issue. However and even though AMPK involvement in LKB1 functions appears until now preponderant, the LKB1 complex is able to activate 13 other kinases belonging to the AMPK family [START_REF] Sapkota | Ionizing radiation induces ataxia telangiectasia mutated kinase (ATM)-mediated phosphorylation of LKB1/STK11 at Thr-366[END_REF] which are poorly studied compare to AMPK ( ). Most of these kinases have been confirmed to be activated by LKB1 in figure 2 physiological conditions since their activity is significantly reduced in Lkb1-deficient models.

Among them, the two brain-specific kinases ( ) also called SAD kinases (SAD-B and A), mainly expressed in brain BRSK-1 and -2 and testis, play a role in neuronal cell polarity. BRSK deficient mice exhibits drastic reduced axon growth as well as mislocalized axons and dentritic markers , . In addition, BRSK-2 has been demonstrated to regulate centrosome duplication. Indeed, BRSK-2 localizes [75 76 ] at centrosome through its association with Tubulin and phosphorylates it at serine 131 which is required to centrosome duplication .

γ- [START_REF] Alvarado-Kristensson | SADB phosphorylation of gamma-tubulin regulates centrosome duplication[END_REF] Thus, BRSK-2 could regulate polarity and neuronal proliferation as AMPK does for epithelial cells (see below). Moreover, like AMPK, BRSK-2 can be also phosphorylated and activated, independently from LKB1, notably by CAMKK .

[78 ]

The LKB1 complex had been described to phosphorylate the four human kinases (microtubule affinity regulating kinase), MARK homologous to the gene product . Even though, MARK kinases appear involved in epithelial cell polarity , , C. elegans Par-1 [START_REF] Lizcano | LKB1 is a master kinase that activates 13 kinases of the AMPK subfamily, including MARK/PAR-1[END_REF] [ 79 80 ] the relationship between LKB1 and MARK kinases remain unclear in mammals. Indeed, LKB1 deficiency either in primary fibroblasts or in cardiac muscle cells does not lead to significant modifications of MARK3 phosphorylation and activation . This could be explained [START_REF] Lizcano | LKB1 is a master kinase that activates 13 kinases of the AMPK subfamily, including MARK/PAR-1[END_REF] by the fact that TAO1 kinase (Thousand and one amino acid protein kinase) is able to phosphorylate the critical T-loop threonine and thereby activate MARK kinases . Nevertheless, forced expression of active LKB1 triggers MARK2 function by suppressing tubulin [START_REF] Timm | MARKK, a Ste20-like kinase, activates the polarity-inducing kinase MARK/PAR-1[END_REF] polymerization through Tau phosphorylation and affecting cell migration and polarity . Thus, in physiological conditions, MARK [START_REF] Kojima | Suppression of tubulin polymerization by the LKB1-microtubule-associated protein/microtubule affinity-regulating kinase signaling[END_REF] kinase activation by LKB1 could be restricted to some isoforms.

Like for other related AMPK kinases, activation of and kinases is due to their phosphorylation on the T-loop threonine by QSK SIK LKB1. For these kinases, this motif is a docking site for the 14-3-3 family members which participates to their full activation and concomitantly affects their subcellular localization . Indeed, an active form of QSK localizes in punctate structures of unknow nature

[41 ]
within the cytosol through 14-3-3 binding . While QSK functions are not well established, QSK knockdown in triggers

[41 ]
Drosophila mitotic defects like spindle and chromosome alignement abnormalities indicating that this protein might be involved also in cell proliferation . In addition to the indication that both AMPK and LKB1 deficiencies lead to similar defects , , it is detachment, called ano kis . This study shows that following cell detachment, SIK is activated in a LKB1 dependent pathway ï [START_REF] Cheng | SIK1 couples LKB1 to p53-dependent anoikis and suppresses metastasis[END_REF] triggering p53 phosphorylation and apoptosis through proapoptotic Bax and Puma induced expression. Thus, loss of LKB1, SIK1 or p53 interferes in ano kis induction so favouring metastasis development .

, also called QIK, appears to be highly expressed in ï [START_REF] Cheng | SIK1 couples LKB1 to p53-dependent anoikis and suppresses metastasis[END_REF] SIK2 adipose tissues and is able to phosphorylate insulin receptor substrate 1 (IRS1) at the same site than AMPK . Similarly, SIK2 [START_REF] Screaton | The CREB coactivator TORC2 functions as a calcium-and cAMP-sensitive coincidence detector[END_REF] phosphorylates the CREB coactivator TORC2 at the same site than AMPK, leading to TORC2 capture by 14-3-3 proteins. Result of this binding is TORC2 retention in the cytosol preventing its transcriptional induction of gluconeogenic genes expression .

[85 ]

About

, also referenced as ARK5, few things are known, except that mice deficient for this kinase exhibit an embryonic NUAK2 lethality after 16.5 days of development due to an exencephaly which results in a brain protrusion outside the skull . Moreover, [START_REF] Tsuchihara | Susceptibility of Snark-deficient mice to azoxymethane-induced colorectal tumorigenesis and the formation of aberrant crypt foci[END_REF] NUAK2 has been described to be activated through TNF , leading to inhibition of the actomyosin contractile network through α phosphorylation of threonine 696 and 853 of MYPT1 regulatory subunit of the myosin light chain phosphatase. Thereby, this inhibition induces increased phosphorylation of myosin light chain and cellular contractility . Unfortunately, this study does not link this TNF [START_REF] Yamamoto | Identification of a novel substrate for TNFalpha-induced kinase NUAK2[END_REF] αdependent NUAK2 activity with the LKB1 complex. Nevertheless, activated by LKB1 complex has been found recently to lead NUAK1 to similar effect on MYPT1

. This suggests that, through NUAK kinases, LKB1 could play an important role in cell contractility [START_REF] Zagorska | New roles for the LKB1-NUAK pathway in controlling myosin phosphatase complexes and cell adhesion[END_REF] regulation independently of Rho GTPAses family effector such as Rho kinase . Besides, NUAK1 has been reported to be increased in [START_REF] Yamamoto | Identification of a novel substrate for TNFalpha-induced kinase NUAK2[END_REF] colorectal tumours where it promotes invasiveness . NUAK1 is also involved in cell senescence and apoptosis inhibition .

[89 ] [START_REF] Humbert | Regulation of ploidy and senescence by the AMPK-related kinase NUAK1[END_REF] [91 ]

Beyond an already established connection between cell contractility and apoptosis , contractility driven by NUAK1 could have a [START_REF] Hebert | Rho-ROCK-dependent ezrin-radixin-moesin phosphorylation regulates Fas-mediated apoptosis in Jurkat cells[END_REF] major role in tumorogenesis progress in addition to these well established roles in invasion and metastatic processes.

Finally, , the less studied member of AMPK related kinases, is mainly expressed in testis and might be involved in SNRK spermatogenesis or spermatozoid motility . This idea is coherent with the fact that male mice which express low levels of LKB1 are [START_REF] Jaleel | Identification of the sucrose non-fermenting related kinase SNRK, as a novel LKB1 substrate[END_REF] sterile .

[94 ]

LKB1 and cell polarity

From the seminal studies in , that characterized PAR4, the LKB1 homologue and its involvement in asymmetric cell C.elegans division, implication of this kinase in cell polarity and cell division had been suspected

. LKB1 contribution in cell cycle regulation was [START_REF] Kemphues | Identification of genes required for cytoplasmic localization in early C. elegans embryos[END_REF] confirmed by ectopic expression in various cellular models leading to induction of p21waf Cdk inhibitor in a p53 dependant manner , [START_REF] Tiainen | Growth arrest by the LKB1 tumor suppressor: induction of p21(WAF1/CIP1)[END_REF] . Beyond this effect on cell cycle regulation, the LKB1 complex could be directly involved in cell division process, as its knock down 32 ] in S2 cells affects chromosome alignment, spindle formation and centrosome number . These effects could be the result of Drosophila [START_REF] Bettencourt-Dias | Genome-wide survey of protein kinases required for cell cycle progression[END_REF] loss of either QSK, MARK-3 or AMPK activities which are individually important in cell division . Nevertheless, those findings [START_REF] Bettencourt-Dias | Genome-wide survey of protein kinases required for cell cycle progression[END_REF] remain fragmentary and need to be confirmed in mammalian models.

In addition to PAR4, Kemphues and colleagues had characterized other genes involved as partitioning defective mutants that include mammalian PAR3 and PAR6 co-homologues. These proteins are well established cell polarity proteins and act together with PAR1 and PAR5, homologues of MARK kinases and 14-3-3 respectively. Thus, the role of LKB1 in cell polarity establishment was initially investigated and observed in oocyte polarization . In mammals, LKB1 complex contribution to apico-basal cell polarity Drosophila [START_REF] Martin | A role for Drosophila LKB1 in anterior-posterior axis formation and epithelial polarity[END_REF] was revealed by co-expressing STRAD and LKB1 in colorectal cancer cell lines . Indeed, in those cells, forced LKB1 activity α [START_REF] Baas | Complete polarization of single intestinal epithelial cells upon activation of LKB1 by STRAD[END_REF] induced leads to formation of brush border-like structures and an actin cap resembling to the apical domain of polarized epithelial cells [START_REF] Baas | Complete polarization of single intestinal epithelial cells upon activation of LKB1 by STRAD[END_REF] . Moreover, this apical-like domain was concomitantly associated with directed endosomal transport of basolateral markers, such as the ] transferin receptor, excluded from the apical structure . However, LKB1 involvement in epithelial cell polarity is less obvious than [START_REF] Baas | Complete polarization of single intestinal epithelial cells upon activation of LKB1 by STRAD[END_REF] usually considered. Indeed, several reports describe that LKB1 knock-down in MDCK , or CaCo2 cells derived respectively

[53 61 ] [53 ]
from canine kidney or human colon, only exhibits a slight delay in cell polarity. Moreover, during embryonic development of LKB1 deficient mice, no cell polarity defect is observed until 7 days of development while important processes of cell polarization have already taken place before this stage . Studies in can potentially explain this discrepancy. Indeed, it has been shown that AMPK or [START_REF] Ylikorkala | Vascular abnormalities and deregulation of VEGF in Lkb1-deficient mice[END_REF] Drosophila LKB1 loss leads to defects in epithelial cell polarity only under energetic stress . These results suggest that LKB1-AMPK Drosophila [START_REF] Mirouse | LKB1 and AMPK maintain epithelial cell polarity under energetic stress[END_REF] is required to allow cell adaptation to unfavourable growth conditions, probably because epithelial cell polarity maintenance is an energetic consuming process. It would be interesting to define whether, in mammalian cells, a similar mechanism is conserved. Indeed, during epithelial derived tumours formation loss of polarity could be due not only to loss of structural protein functions but result also to limiting nutriment availability leading to an energetic stress that tumour cells have difficulties to bypass. Thus, observations that E-cadherin is partly involved in LKB1 dependent AMPK activation and often lost during tumour progression , could reflect another way for it [53 ] [95 ] to affect cell polarity in addition to its well established structural functions . This eventuality is strengthened by the fact that AMPK [START_REF] Cavallaro | Cell adhesion and signalling by cadherins and Ig-CAMs in cancer[END_REF] activity, involved in tight junction formation , , is increased in epithelial cells through calcium switch , considered in part to [97 98 ] [97 98 ] act on E-cadherin and adherens junctions formation. Main role of the tight junctions is to seal the apical membranes of adjacent epithelial cells precluding the free diffusion of solutes toward the basolateral membranes . However, they also restrict the lateral diffusion of [START_REF] Caudron | Septins and the lateral compartmentalization of eukaryotic membranes[END_REF] some lipids and plasma membrane-associated proteins between apical and basolateral domains, which is thought to ensure that the two membrane domains do not mix . How AMPK activity controls these events is unclear. It seems that its inhibitory properties on the [START_REF] Caudron | Septins and the lateral compartmentalization of eukaryotic membranes[END_REF] mTOR pathway may play a role, since the effect of dominant negative AMPK construct on tight junctions formation can be partially rescued by mTOR inhibitors such as Rapamycin . This is consistent with the fact that mTOR activity is opposed to apico-basal [START_REF] Zhang | AMP-activated protein kinase regulates the assembly of epithelial tight junctions[END_REF] polarity by contributing to the epithelio-mesenchymal transition . All of this underlines that outside of energetic stress conditions where LKB1-AMPK appears to be dispensable in epithelial cell polarity, the mTOR pathway activity should be maintained low by an alternative mechanism. Moreover, MARK kinase activities, well known to be required in epithelial cell polarity , would be also [79 80 ] activated through a LKB1 independent manner potentially by TAO1 kinase as mentioned above .

[81 ]

Conclusions

Currently, there is no doubt that the LKB1 complex is at the centre of an important signalling node affecting numerous cellular processes whose deregulation contributes to pathologies such as Peutz-Jeghers syndrome and cancer. Although studies by many labs have focused on the role of LKB1 complex in the regulation of the AMPK-mTOR pathway, the 13 AMPK related LKB1 effectors remain currently poorly investigated. Thus, studies using LKB1 knock down models and relying consequence almost exclusively on the AMPK activity should be interpreted with caution since contribution of other kinases of this family might participate to consequences observed.

Besides and although downstream signalling events of this complex begin to be defined, upstream signals that regulate LKB1 functions remain almost completely uncharacterized. Expressions of MO25 and STRAD paralogues as well as splice variants have probably some importance in specific LKB1 functions. The emerging notion that the complex is constitutively active in cells led to think that the subcellular localizations of LKB1 are probably a major feature to specify the functions of the kinase. How these localizations are achieved and which effectors and substrates are targeted are pending questions that should be addressed in a near future to understand all aspects of LKB1 complex involvement and disorders in physiopathology.

85 ]

 85 tempting to speculate that those three kinases may cooperate in a same pathway during cell division, potentially linked to Drosophila energetic stress. salt-induced kinase 1) has been described to play a role in steroidogenesis and TGF signalling, for this latter by inhibiting SIK1 ( β type I TGF receptor kinase . More recently, SIK1 was also found to have a crucial role in p53 dependent apoptosis induced by cell β[START_REF] Takemori | Salt-inducible kinase-mediated regulation of steroidogenesis at the early stage of ACTH-stimulation[END_REF]

′Fig. 1 (

 1 Fig. 1 (A) Scheme of human LKB1 spliced variants, STRAD and Mo25 paralogues. Tilted and underline residues considered as autophosphorylation sites. Phosphorylated residues targeted by mentioned kinases (arrow). NLS for nuclear localization sequence. CAAX for the site of farnesylation. (B) Schematic representation of active and functional complex.

Fig. 2

 2 Fig. 2 Domain organization of AMPK-related kinases activated by LKB1 complexes. Kinase domain in light gray, ubiquitin associated domain in checked motif. Usual and alternative names, molecular weight and main identified functions, respectively from left to right.

Acknowledgements:

We apologize to researchers whose investigations were not included in this review, due to space limitations. This work is supported by INSERM and La Ligue Nationale Contre le Cancer (Label Ligue 2010 JPB).


José V Manjón
email: jmanjon@fis.upv.es

Pierrick Coupé
email: pierrick.coupe@gmail.com

Antonio Buades
email: toni.buades@uib.es

D Louis Collins

Montserrat Robles
email: mrobles@fis.upv.es

New Methods for MRI Denoising based on Sparseness and Self-Similarity

Keywords: MRI, denoising, sparseness, nonlocal means VBM: Voxel-based morphometry PCA: Principal component analysis DCT: Discrete cosine transform NLM: Nonlocal means ODCT3D: Oracle-based 3D discrete cosine transform filter RI-NLM3D: Rotationally invariant nonlocal means filter PRI-NLM3D: Prefiltered rotationally invariant nonlocal means filter WSM: Wavelet sub-band mixing RMSE: Root mean squared error SSIM: Structural similarity index

This paper proposes two new methods for the three-dimensional denoising of magnetic resonance images that exploit the sparseness and self-similarity properties of the images. The proposed methods are based on a threedimensional moving-window discrete cosine transform hard thresholding and a three-dimensional rotationally invariant version of the well-known nonlocal means filter. The proposed approaches were compared with related state-ofthe-art methods and produced very competitive results. Both methods run in less than a minute, making them usable in most clinical and research settings.

Introduction

The denoising of magnetic resonance (MR) images is an important open problem that has been discussed frequently in the recent literature because of its relevance to many clinical and research uses. Denoising is used as a preprocessing step in many image processing and analysis tasks such as registration or segmentation to reduce the random noise arising from the acquisition process.

One approach that has been applied extensively in MRI preprocessing is the Gaussian filter [START_REF] Ashburner | Voxel-based morphometry-the methods[END_REF]. This method, although capable of reducing some image noise (especially in homogeneous areas), also removes high-frequency signal components, thereby blurring edges in the images.

Therefore, this filter has been commonly used for regularization purposes, such as in voxel-based morphometry (VBM) [START_REF] Ashburner | Voxel-based morphometry-the methods[END_REF], to reduce anatomical inconsistencies.

A large number of edge-preserving methods have been proposed to overcome the above-mentioned blurring effects. For example, anisotropic diffusion filters [START_REF] Gerig | Nonlinear anisotropic filtering of MRI data[END_REF] are able to remove noise using gradient information while respecting important image structures. Recently, [START_REF] Krissian | Noise-driven anisotropic diffusion filtering of MRI[END_REF] proposed a new anisotropic diffusion filter based on a linear minimum mean squared error estimation and partial difference equations for Rician noise removal that has achieved state-of-the-art results.

Wavelet-based filters have also been applied successfully to the denoising of MR images [START_REF] Pizurica | A versatile wavelet domain noise filtration technique for medical imaging[END_REF]. Such filters are rooted in the processing of images in a transformed domain. Other transforms that have been applied to denoise images include principal component analysis (PCA) [START_REF] Muresan | Adaptive principal components and image denoising[END_REF] and the discrete cosine transform (DCT) [START_REF] Yaroslavsky | Transform domain image restoration methods: review, comparison and interpretation[END_REF]. Many transform domain filters derive from variations of the transformthreshold-inverse transform principle.

In accordance with this principle, local transform approaches (i.e., sliding window with or without overlapping) have obtained very good results in recent years [START_REF] Guleryuz | Weighted overcomplete denoising[END_REF][START_REF] Guleryuz | Weighted averaging for denoising with overcomplete dictionaries[END_REF][START_REF] Yaroslavsky | Transform domain image restoration methods: review, comparison and interpretation[END_REF]. In [START_REF] Guleryuz | Weighted averaging for denoising with overcomplete dictionaries[END_REF] method for Gaussian noise reduction, the image noise is removed using overcomplete linear transforms and thresholding. In practice, Guleryuz applied a classical sliding-window DCT thresholding as in [START_REF] Yaroslavsky | Transform domain image restoration methods: review, comparison and interpretation[END_REF], but overlapping estimations were adaptively combined to reduce the Gibbs effects.

Other recently proposed approaches use learned image patch dictionaries [START_REF] Aharon | K-SVD: an algorithm for designing overcomplete dictionaries for sparse representation[END_REF][START_REF] Elad | Image denoising via learned dictionaries and sparse representation[END_REF][START_REF] Mairal | Learning multiscale sparse representations for image and video restoration[END_REF] instead of DCT bases to perform the denoising. All these approaches stem from the fact that an image can be represented as the linear combination of a set of image bases with very few non-null coefficients. This property, known as sparseness, is the core of the JPEG and JPEG2000 compression standards.

Finally, the nonlocal means (NLM) filter, a new method introduced by [START_REF] Buades | A non-local algorithm for image denoising[END_REF], has emerged as a very simple and effective way to reduce noise while minimally affecting the original structures of the image. This method is based on the natural redundancy of patterns within the images. Recently improved, the NLM filter has been applied to the denoising of MR images and demonstrated better results than those of previous related methods (Coupé et al., 2008a(Coupé et al., , 2008b;;[START_REF] Manjón | Adaptive non-local means denoising of MR images with spatially varying noise levels[END_REF]Wiest-Daesllé et al., 2008).

In this paper, we present two new approaches to the three-dimensional (3D) denoising of MR images. The first is an extension of the original method proposed by [START_REF] Guleryuz | Weighted averaging for denoising with overcomplete dictionaries[END_REF]. Based on local 3D DCT hard thresholding, our proposed method has been adapted to deal with Rician noise (the typical type of noise present in magnitude MR images) using a pseudo-oracle principle. The second proposed method is a new rotationally invariant 3D version of the Rician-adapted NLM filter that uses a prefiltered image obtained by DCT denoising.

Methods

Typically, the observed noisy MR image y is considered the linear combination of a noise-free image x and a white noise realization n from the measurement process:

n x y   (1)
Therefore, the goal of any denoising algorithm is to find a good estimate x ˆ, given y. This section describes two recent methods for image denoising based on efficient concepts of how to solve this problem:



Denoising using sparseness: The first method under study is the local DCT denoising method proposed by [START_REF] Guleryuz | Weighted averaging for denoising with overcomplete dictionaries[END_REF] that uses the sparseness of the image (i.e., the ability of the image to be represented by a small number of base functions, e.g., DCT bases).



Denoising using self-similarity:

The second method under investigation is the well-known NLM filter proposed by [START_REF] Buades | A non-local algorithm for image denoising[END_REF] that uses the pattern redundancy present in the image. After a brief presentation of the original NLM filter, a new rotationally invariant version of this filter is proposed.

In order to take advantage of both properties, sparseness and self-similarity, a new efficient collaborative method is proposed. The proposed technique is based on a rotationally invariant version of the NLM filter and uses an image prefiltered with the extended DCT-based method to compute patch similarities.

Finally, the adaptation of the proposed methods to allow them to handle Rician noise (i.e., the type of noise present in magnitude MR images) is described.

Denoising using sparseness

Guleryuz's method (DCT3D)

In [START_REF] Guleryuz | Weighted averaging for denoising with overcomplete dictionaries[END_REF] method, the estimate of the noise-free image x ˆis obtained using a hard thresholding method. An overcomplete set of 3D block DCTs (444 block size) is used (Guleryuz used 88 blocks in his 2D method, which yields the same number of coefficients). The local denoised estimate at block j,

j
x ˆ, is obtained by applying a hard thresholding rule:

) ( j j y c   ) , ( ˆ j j c T c  ) ( ˆ1 j j c H x   (2)
where H is a 3D DCT, c j are the transform coefficients of the j block, and T is the hard thresholding operator with threshold . Finally, all local estimates ) ( ˆi x j are combined from all overlapping j blocks at position i using the following weighted average rule:

     V j j V j j j i x i x 1 1 ) ( ) ( ˆ  0 1 1 j j c    (3)
where V is the number of overlapping blocks contributing to  ˆ x ( i) and  j is the weight of each block j, which is proportional to the inverse of the j cˆ L0 norm (number of nonzero coefficients of block j after the thresholding operation). This approach gives more weight to estimations with more null values after thresholding. Thus, homogeneous blocks tend to be favored in comparison with blocks containing edges and thus suffering from Gibbs effects. The only parameter in this method is the threshold value .

Oracle-based filters

Oracle-based filters assume that, if the null coefficients of the original noise-free image are known, they can be applied to improve the denoising. In fact, the thresholding method that minimizes the mean squared error with the original image is given by:

) ( j j x H p         otherwise i c i p if i c j j j ) ( ) ) ( ( 0 ) ( ˆ (4)
which compares the magnitude of the noise-free DCT coefficients with the standard deviation of the noise  (see [START_REF] Mallat | A wavelet tour of signal processing[END_REF] for more details).

Unfortunately, in practice, no noise-free image x is available; hence, this approach is unfeasible.

Proposed Oracle-based DCT filter (ODCT3D)

If we relax the Oracle condition, we can use a prefiltered image pre x ˆ with, for example, the described DCT3D method as an approximation of the noise-free image. With this approach, the c j block coefficients of the noisy image y can be better thresholded (applying the threshold operation described in Equation 4) using the corresponding block coefficients p j of the prefiltered image. The final image reconstruction is performed using the same approach described above.

We will refer to this method as Oracle DCT3D (ODCT3D). As we will show later, this modification not only allows better results to be obtained in terms of quantitative measures, but also reduces the remaining Gibbs artifacts of the original DCT3D method.

The only parameter of this method is the threshold  from the prefiltering step.

In all of our experiments, we used a  value of 2.7  (where  is the standard deviation of the noise), which is common in DCT and wavelet thresholding methods [START_REF] Mallat | A wavelet tour of signal processing[END_REF].

Denoising using self-similarity

The original nonlocal means denoising filter

For its use of the self-similarity properties of an image, the NLM filter has recently become known as a simple but effective way of removing noise from an image while minimally affecting its original structure.

Originally proposed by [START_REF] Buades | A non-local algorithm for image denoising[END_REF], the NLM filter takes advantage of the high level of pattern redundancy in an image, achieving high-quality image denoising by averaging similar realizations of the noisy signals. Basically, this filter reduces the noise in an image by averaging voxels that originally had the same intensity in the noise-free image. To this end, [START_REF] Buades | A non-local algorithm for image denoising[END_REF] suggested that voxels with similar neighborhoods (small 3D patches in our volumetric case) tend to have similar original values, yielding the following filter expression:

       j j j i w i y j i w i x ) , ( ) ( ) , ( ) ( ˆ 2 2 2 ) , ( h N N j i e j i w    (5)
where  represents the search volume, the weights w(i,j) represent the similarity between any two 3D patches N i and N j centered around pixels i and j, and h 2 is a parameter controlling the strength of the filter. So defined, this filter is able to identify N j patches as similar to a given patch N i , such that they share the same structure and orientation. This means that the weights defined in Equation 5 are not rotationally invariant. Therefore, similar patches with similar structure but different orientations to the reference patch will have a small influence in the average.

Rotationally invariant nonlocal means filter (RI-NLM3D)

To avoid this problem, the orientation of each patch can be estimated and corrected before computing the L2 norm [START_REF] Kleinschmidt | Nonlocal texture filtering with efficient tree structures and invariant patch similarity measures[END_REF][START_REF] Zimmer | A rotationally invariant block matching strategy improving image denoising with non-local means[END_REF]. However, this solution is very time consuming, as it increases the already high computational burden of the NLM method, and the reorientation operation introduces interpolation artifacts and correlates the noise, making the denoising process more difficult.

Instead of correcting the orientation of the patches, a rotationally invariant measure describing the patch information can be used [START_REF] Lou | Nonlocal similarity image filtering[END_REF][START_REF] Zexuan | A moment-based nonlocal-means algorithm for image denoising[END_REF], for instance, rotationally invariant image descriptors such as Hu or Zernike moments [START_REF] Hu | Visual pattern recognition by moment invariants[END_REF][START_REF] Teague | Image analysis via the general theory of moments[END_REF]. However, such descriptors are sensitive to noise, and their ability to represent patch information decreases rapidly as noise power increases.

Therefore, we present a simple but effective similarity measure that is rotationally invariant and based on voxel intensity and the corresponding local patch mean (using a 333 voxel Gaussian kernel):

             2 2 2 2 ) ( 3 )) ( ) ( ( 2 1 ) , ( h j y i y Nj i N e j i w   (6)
where µ Ni and µ Ni are the mean values of patches N i and N j around voxels i and j, and h is related to the standard deviation of the noise. Note that the distance between the patch means is multiplied by three to account for this distance being approximately three times smaller than the average voxel distance because of the application of Gaussian smoothing. Our experiments showed that giving this distance three times the weight produces the best results.

Because it uses the value of the central voxel of the patch, this similarity measure is point specific rather than region specific, which allows a better description of the characteristics of the voxel. On the other hand, using the mean value of the region around the voxel would make the similarity measure more robust against the noise while minimizing the aggregation artifacts that could produce cartoon-like images.

In analyzing Equation 6, it becomes clear that this measure exhibits the same noise sensitivity as Hu or Zernike moments because it uses the intensity of the voxel of the similarity measure, which is not a robust descriptor of the real intensity of the voxel (the main benefit of this approach is its low computational complexity). For this reason, applying this filter to medium and high noise levels produces poor results. However, we will show that, under certain circumstances, this method can obtain very good denoising results when combined with a prefiltering step. We will refer to this method as the RI-NLM3D.

A new collaborative method based on sparseness and self-similarity

We propose to combine the two methods described above to obtain better results than could be obtained by either method separately. Specifically, we propose an Oracle-based RI-NLM3D method where the similarities between voxels and patches are computed from the already denoised image using the ODCT3D method and then applied to denoise the original noisy image:

       j j j i i y j i i x ) , ( ) ( ) , ( ) ( ˆ                 2 2 2 2 ) ( 3 )) ( ) ( ( 2 1 ) , ( h j x i x o Nj o i N o o e j i    (7)
where the new weight  is computed using the ODCT3D denoised image x o and applied to denoise the noisy image y. [START_REF] Kervrann | Bayesian non-local means filter, image redundancy and adaptive dictionaries for noise removal[END_REF] proposed a similar plug-in estimator for 2D NLM filtering.

In their approach, they first denoised the image with a Bayesian NLM filter. The similarity between denoised patches was then computed and finally applied to noisy voxels. However, this method has the disadvantage of being computationally intensive and only based on pattern redundancy without taking into consideration the sparseness properties of the image.

By contrast, the proposed approach presents an interesting advantage in that the RI-NLM3D method only works well for low noise levels, a condition that is met after the ODCT3D method has been applied. Figure 1 shows a scheme of the proposed method. Finally, a preselection strategy was also applied to avoid useless computations (Coupé et al., 2008a;[START_REF] Kervrann | Bayesian non-local means filter, image redundancy and adaptive dictionaries for noise removal[END_REF]. We used a statistically driven rule based on the distance between the patch means of the prefiltered volume:

                      . 0 ) , ( 2 2 2 2 ) ( 3 )) ( ) ( ( 2 1 otherwise h if e j i o j o i h j x i x o Nj o i N o o      (8)
In this way, patches with an intensity difference greater than h do not contribute to the denoising process. We will refer to this proposed method as the prefiltered rotationally invariant NLM3D (PRI-NLM3D).

Because the proposed similarity measure is applied to the already denoised image using the ODCT3D method, the h parameter controlling the strength of the filter was set to 0.4  instead of  (this value was found experimentally to be the best option through an exhaustive search of many image types and noise conditions). As in Coupé et al. (2008a), the radius of the search volume  was set to 5 voxels, which resulted in a 3D search volume of 111111 voxels.

Experimentally, this volume was shown to be a good compromise between computational complexity and quality of the results. Larger search volumes improved the results only marginally, but increased the computational burden significantly.

Adaptation to Rician noise

Noise in magnitude MR images follows a Rician distribution [START_REF] Nowak | Wavelet-based Rician noise removal for magnetic resonance imaging[END_REF].

Consequently, the weighted average produced by applying the PRI-NLM3D method will be biased because of the asymmetry of the Rician distribution.

To avoid such bias, we adopted a method of bias correction used in early MRI denoising [START_REF] Manjón | MRI denoising using non-local means[END_REF][START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio MRI: applications to DR-MRI[END_REF]. The unbiased PRI-NLM3D estimator is defined as follows:

                              0 , 2 ) , ( ) ( ) , ( max ) ( ˆ 2 2    i i j j j i i y j i i x (9)
We also adapted the ODCT3D method to deal with Rician noise. In this case, because of the effect of DCT thresholding, the bias in the squared domain is not constant, but dependent on intensity. It can be estimated theoretically and inverted in the original domain using the properties of the first moment of a Rician distribution, as expressed in Equation 10:

2 2 1 2 2 0 2 2 j 4 ) ( 2 ) ( 4 ) ( 2 ) ( 1 2 ) ( exp 2 (i) x ˆ                                                      i I i i I i i j j j j j       (10) with   ) ( ) ( i x i j j  where ) ( ˆi x j and ) ( i x j
are the biased and unbiased estimates due to the Rician bias at position i from block j, I 0 and I 1 are the modified Bessel functions of order zero and one, respectively, and  is the standard deviation of the noise. Then, the final unbiased estimate is obtained by mapping the biased estimates directly onto the unbiased ones by using a precomputed look-up table. Finally, the different estimations from the different blocks are averaged to provide the final unbiased estimate in a manner similar to Equation 3:

     V j j V j j j i x i x 1 1 )) ( ( ) ( ˆ    (11)
where operator  is a mapping function that provides the unbiased estimate for a given biased estimate by using the precomputed equivalences from Equation 10.

Experiments and results

To evaluate the described methods, we used the well-known BrainWeb 3D MRI phantoms [START_REF] Collins | Design and construction of a realistic digital brain phantom[END_REF][START_REF] Kwan | MRI simulation-based evaluation of image-processing and classification methods[END_REF], T1-weighted (T1w), PDw, and T2w volumes of 181217181 voxels (voxel resolution = 1 mm 3 ), which were corrupted with different levels of Rician noise (1% to 15% of maximum intensity). Rician noise was generated by adding Gaussian noise to real and imaginary parts and then computing the magnitude image.

Two quality measures were used to evaluate the results. The first was the root mean squared error (RMSE) metric, which is a distance measure commonly used in image processing. The second was the structural similarity index (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF], which is a measure more consistent with the human visual system:

) )( ( ) 2 )( 2 ( ) , ( 2 2 2 1 2 2 2 c c c y x SSIM y x y x xy y x              , ( 12
)
where µ x and µ y are the mean value of images x and y,  x and  y are the standard deviation of images x and y,  xy is the covariance of x and y, c 1 = (k 1 L) 2 , and c 2 = (k 2 L) 2 (where L is the dynamic range, k 1 = 0.01 and k 2 = 0.03). As suggested by [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF], the SSIM was locally estimated using a Gaussian kernel of 333 voxels. Finally, the mean value of all the local estimations was used as a quality metric. For the sake of clarity, both measures were estimated only in the region of interest (head tissues) obtained by removing the background (i.e., the label 0 of the discrete model in BrainWeb).

DCT3D vs. ODCT3D

To evaluate the proposed ODCT3D method, we compared it with the 3D version of Gurelyuz's ( 2007) method (DCT3D) using the T1w BrainWeb phantom. Figure 2 shows the RMSE and SSIM measures for different levels of Rician noise and indicates that the proposed method outperforms the DCT3D method at medium and high noise levels. Furthermore, as can be seen in Fig. 3, the proposed ODCT3D method reduces the Gibbs effects that remain after the DCT3D method is applied. methods for a Rician noise level of 15%. The DCT3D filtered image shows some Gibbs effects, mainly observable in homogeneous areas, while the image produced with the proposed ODCT3D method no longer shows such artifacts.

Comparison of methods

We compared our proposed Rician-adapted versions of the ODCT3D and PRI-NLM3D methods with some recently proposed related methods used for MRI denoising of Rician distributed noise: the optimized blockwise nonlocal means filter (Coupé et al., 2008a), the wavelet sub-band coefficient mixing method (WSM) (Coupé et al., 2008b), and the ORNRAD filter [START_REF] Krissian | Noise-driven anisotropic diffusion filtering of MRI[END_REF].

Figure 4 shows a comparison of the RMSE and SSIM measures obtained using the methods under investigation with different image types and noise levels.

Figures 5,6, and 7 are provided for a visual comparison of the results. All the methods were run with the optimal parameters suggested by the authors (see Table 1). With the exception of the ORNRAD filter, all the methods were implemented using multithreaded C MEX files in MATLAB 2009a (The MathWorks, Inc.) in Windows XP 64-bit Edition (Pentium quad-core 2.4 GHz with 8 GB of RAM). The ORNRAD filter was run (using AMILab 2.0.4 on an Intel Xeon X5650 with 12 cores and 12 GB of RAM running Linux Fedora 13) by its authors, who supplied results for the same datasets.

To facilitate the reproducibility of the presented experiments, the Matlab code for our the experiments will be available on our webpage: http://personales.upv.es/jmanjon/denoising/prinlm.html.

Table 1. Method parameters (v is the radius of the search volume, f is the radius of the 3D patches, h is the strength of the filter, and  is the standard deviation of the noise).

Method Parameters

Blockwise As can be seen in Fig. 4, the proposed ODCT3D and PRI-NLM3D methods showed improvements over the other filters under consideration in most of the cases. The PRI-NLM3D method was the best method in almost all cases. Only the ODCT3D method surpassed the PRI-NLM3D method at high noise levels on the PDw data, probably because of the low contrast of the images. The ORNRAD method performed very well on the T1w images, producing similar results to those of the proposed ODCT3D method, although its performance decreased noticeably on the PDw and T2w images in terms of the RMSE.

NLM v = 5, f = 1, h =  WSM v = 3, f 1 = 1, f 2 = 2, h =  ODCT3D  = 2.7  PRI-NLM3D v = 5, h = 0.4  ORNRAD v = 1, T = 2, dt = 1/6,  1 = 0.7,  2 = 1
Furthermore, as can be observed in Figs. 5, 6, and 7, the PRI-NLM3D produced visually more pleasant results than the other methods, showing fewer oscillations in homogeneous areas and better defined edges.

In terms of filtering time, the ORNRAD method was the slowest, taking 12 min on average (note that this was using a single-threaded implementation, whereas the other methods were running on four cores), followed by the WSM method at 110 s, the blockwise NLM3D at 81 s, the PRI-NLM3D method at 60 s, and finally the ODCT3D method at only 10 s.

Evaluation on real clinical data

To evaluate the consistency of the PRI-NLM3D method on clinical data, two The Rician noise level was estimated to be around 1% of maximum intensity.

Figure 9 allows for a visual comparison of the results produced using the ODCT3D, ORNRAD, and PRI-NLM3D methods. The PRI-NLM3D method removed the noise successfully while preserving fine details of the images, whereas the ODCT3D and especially the ORNRAD method slightly oversmoothed some details. The processing time for this dataset was 9 s with the ODCT3D method, 52 s with the PRI-NLM3D method, and 13 min with the ORNRAD filter.

Discussion

We presented two new methods for MRI denoising (ODCT3D and PRI-NLM3D) that take advantage of two intrinsic properties of MR images: sparseness and self-similarity. The proposed PRI-NLM3D method was compared with state-ofthe-art methods in MRI denoising and produced the best results out of all the methods under consideration.

The enhanced accuracy of the proposed methods can be understood by taking two main facts into consideration. First, the ODCT3D method benefits from the high compressibility (sparseness) of MRI data, allowing a very efficient noise reduction (the ODCT3D method obtained very good results in comparison with previously proposed NLM-based methods). Second, the PRI-NLM3D method indirectly takes benefit from the sparseness properties of the images by using the prefiltered ODCT3D data as prior. This method showed a very good performance as a result of using the low noise image prior and the increased number of redundant patterns present when using a rotationally invariant similarity measure (compared with the non-rotationally invariant NLM voxelwise version). It is worth noting that, although other nonlocal rotationally invariant methods have been previously proposed in the literature, they were 2D implementations rather than 3D like the one proposed in this paper.

Both the ODCT3D and PRI-NLM3D methods showed an improved performance over previous state-of-the-art MRI denoising methods. The main differences between these two methods exist in the small details: The ODCT3D method tends to slightly oversmooth edges and some fine details (compared with the PRI-NLM3D), probably as a result of its nature as a block processor, whereas the PRI-NLM3D method seems to retain more detail in the denoised image because of its voxelwise processing.

From a practical point of view, the proposed methods run in a reasonable time (less than 1 min) for most research and clinical settings. Moreover, GPU-based implementations can be used to further reduce the processing time. Recently, a GPU-based NLM filter implementation showed a drastic reduction in the computational burden [START_REF] Huang | Non-local means denoising algorithm accelerated by GPU[END_REF], which could allow for an almost real-time version of the proposed filters.

Finally, it is worth noting that the application of the proposed method can clearly benefit not only visual diagnostics, but also quantitative methodologies that rely on good quality data, such as MRI brain tissue segmentation or MR diffusion tensor image processing.

Fig. 1 .

 1 Fig. 1. Scheme of the proposed method.

Fig. 2 .

 2 Fig. 2. Left: Root mean squared error (RMSE) of the ODCT3D and DCT3D methods for different noise levels. Right: Structural similarity index (SSIM) of the ODCT3D and DCT3D methods for different noise levels.

Fig. 3 .

 3 Fig. 3. Example result of the filtering process (axial slice) using the ODCT3D and DCT3D

Fig. 4 .

 4 Fig. 4. Left: RMSE of the compared methods for different image types and noise levels. Right: SSIM of the compared methods for different image types and noise levels.

Fig. 5 .

 5 Fig. 5. Example filtering results for an axial slice of the T1w BrainWeb phantom (Rician noise level of 15%). The third row shows the absolute value of the image residuals for the different methods.

Fig. 6 .

 6 Fig. 6. Example filtering results for an axial slice of the PDw BrainWeb phantom (Rician noise level of 15%). The third row shows the absolute value of the image residuals for the different methods.

Fig. 7 .

 7 Fig. 7. Example filtering results for an axial slice of the T2w BrainWeb phantom (Rician noise level of 15%). The third row shows the absolute value of the image residuals for the different methods.

  datasets were used. The first was an MP-RAGE T1w volumetric sequence acquired on a Siemens 1.5T Vision scanner (Erlangen, Germany) [TR = 9.7 ms, TE = 4 ms, TI = 20 ms, TD = 200 ms, flip angle = 10 º, voxel resolution = 111.25 mm 3 , 256256128 voxels]. The Rician noise level (used as the filtering parameter) was estimated to be around 2% of the maximum intensity using the object-based method proposed by Coupé et al. (2010). The filtering results for this dataset are shown in Fig. 8. The PRI-NLM3D method removed the noise successfully, showing no significant anatomical information on the image residuals. The processing time for this dataset was 42 s. The second dataset was an MP-RAGE T1w volumetric sequence acquired on a Siemens Trio 3T scanner (Erlangen, Germany) [TR = 2300 ms, TE = 2.9 ms, TI = 900 ms, flip angle = 9 º, voxel resolution = 1 mm 3 , 256240176 voxels].

Fig. 8 .

 8 Fig. 8. Example results of the proposed PRI-NLM3D filter on real data (Rician noise level of 2%). The background and part of the face were removed by a defacer program to preserve the anonymity of the subject. From top to bottom: Original noise volume, denoised volume using the proposed method, and the corresponding residuals.

Fig. 9 .

 9 Fig. 9. Example results of filters applied to real data (Rician noise level of 1%). All compared methods worked very well, but the ODCT3D, and especially the ORNRAD method, tended to oversmooth some fine details in the images (see magnification of the white square region in the lower row).

Acknowledgments

We are grateful to Dr. Krissian for providing us with the results of his denoising method for the comparison section. This work has been partially supported by the Spanish Health Institute Carlos III through the RETICS Combiomed, RD07/0067/2001, the Spanish Ministry Science and Innovation through grant TIN2008-04752, and the Canadian Institutes of Health Research Industry Cda (CECR)-Gevas-OE016.