Michaël Roynard
email: michael.roynard@epita.fr

Edwin Carlinet
email: edwin1.carlinet@epita.fr

Thierry 2022 Géraud
email: thierry.geraud@epita.fr

A Modern C++ Point of View of Programming in Image Processing

Keywords: CCS Concepts:, Software and its engineering → Software development techniques;, Computing methodologies → Image processing Image processing, Generic Programming, Modern C++, Software, Performance ACM Reference Michaël Roynard, Edwin Carlinet,

C++ is a multi-paradigm language that enables the programmer to set up e cient image processing algorithms easily. This language strength comes from many aspects. C++ is high-level, so this enables developing powerful abstractions and mixing di erent programming styles to ease the development. At the same time, C++ is low-level and can fully take advantage of the hardware to deliver the best performance. It is also very portable and highly compatible which allows algorithms to be called from high-level, fast-prototyping languages such as Python or Matlab. One fundamental aspects where C++ shines is generic programming. Generic programming makes it possible to develop and reuse bricks of software on objects (images) of di erent natures (types) without performance loss. Nevertheless, conciliating genericity, e ciency, and simplicity at the same time is not trivial. Modern C++ (post-2011) has brought new features that made it simpler and more powerful. In this paper, we focus on some C++20 aspects of generic programming: ranges, views, and concepts, and see how they extend to images to ease the development of generic image algorithms while lowering the computation time.

Introduction

C++ claims to "leave no room for a lower-level language (except assembler)" [START_REF] Stroustrup | Evolving a Language in and for the Real World: C++ 1991-2006[END_REF] which makes it a go-to language when developing high-performance computing (HPC) image processing applications. The language is designed after a zero-overhead abstraction principle that allows us to devise a high-level but e cient solution to image processing problems. Other aspects of C++ are its stability, its portability on a wide range of architectures, and its direct interface with the C language, which makes C++ easily interoperable with highlevel prototyping languages. Therefore, many image processing libraries (and numerical libraries in general) implement performance-sensitive features in C++ (or C/Fortran as in OpenCV [START_REF] Bradski | The OpenCV Library[END_REF], IPP [START_REF] Burylov | Intel Performance Libraries: Multi-Core-Ready Software for Numeric-Intensive Computation[END_REF]) or with a hardware-dedicated language (e.g. CUDA [START_REF] Brill | NVIDIA VisionWorks toolkit[END_REF]), which are exposed through a high-level API to Python, LUA, etc.

Apart from the performance considerations, the problem lies in that each image processing eld comes with its own set of image type to process. The most common image type is an image of RGB or gray-level values, encoded with 8-bits channel, on a regular 2D rectangular domain. That covers 90% of common usages. However, new image types have come with the development of new devices: 3D multi-band images in Medical Imaging, hyperspectral images in Astronomical Imaging, images with complex values in Signal Processing. An image processing library able to handle those images [START_REF] Stepanov | The standard template library[END_REF] where one has to work on a collection of data. For example, let us consider the algorithm maxof(Collection c) that gets the maximal element of a collection (see g. 4). Whether the collection is actually implemented with a linked-list, a contiguous bu er of elements or whatever data structure is irrelevant. The only requirements for this algorithm are: [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems[END_REF] we can iterate through it; (2) the type of the elements is regular (i.e. behaves the same way as a primitive type like int) and forms a monoid with an associative operator "max" and a neutral element "std::numeric_limits<value_t<R»::min()". A monoid is an abstract algebraic structure equipped with an associative binary operation and an identity element. Such structure is studied to achieve high e ciency in several areas, such as distributed computation of small scale algorithms [START_REF] Dean | Identifying Monoids: Exploiting Compositional Structure in Code[END_REF] calculus of object query languages [START_REF] Fegaras | Towards an E ective Calculus for Object Query Languages[END_REF]. There exists other algebraic structures, such as semi-ring dictionaries [START_REF] Shaikhha | Functional collection programming with semi-ring dictionaries[END_REF] that are studied to unify performance optimizations used in both database and linear algebra.

The constraint (1) is abstracted by pairs of iterators in the STL and ranges in C++20, while C++20 introduces concepts to check if a type follows the requirements of the algorithm. The term concept is de ned as follows in [START_REF] Dehnert | Fundamentals of Generic Programming[END_REF]: "a set of axioms satis ed by a data type and a set of operations we can perform on it.". Concepts enable checking the requirements for (2).

The Image Concept

Most image processing algorithms are also generic [START_REF] Levillain | Why and How to Design a Generic and E cient Image Processing Framework: The Case of the Milena Library[END_REF][START_REF] Levillain | Practical Genericity: Writing Image Processing Algorithms both Reusable and E cient[END_REF][START_REF] Ritter | Image Algebra: An Overview[END_REF] by nature. We saw in section 2.1 that concepts emerges from pattern behavior extracted from algorithms. Similarly to g. [START_REF] Bradski | The OpenCV Library[END_REF] a at structuring element (SE) B (we note Į the SE centered in ). The dilation is de ned as Ĝ ( ) = sup{ ( ), ∈ Į }; the generic algorithm is given in g. 5. As one can see, the implementation does not rely on a speci c implementation of the image. It could be a 2D image, a 3D image or even a graph image (the SE would be the adjacency relation of the graph).

The image requirements can be extracted from this algorithm. The image must provide a way to access its domain which must be iterable. The structuring element must act as a function that returns a range of elements having the same type as the domain element (let us call them points of type ). Image needs to provide a way to access its value for a given point (f(x)) with x of type P. Last, as in g. 4, image values (of type ) have to support max and have a neutral element "0". We deduce the -simpli ed-Image concept and the constrained dilation algorithm in g. 6. Actually, the requirements for being an image are quite light. This provides versatility and allows us to pass non-regular "image" objects as inputs such as the image views in section 3.

While C++20 provides all the tools necessary to properly de ne concepts as well as leveraging them when implementing algorithms, it is still necessary to make the inventory of the algorithms families (explained in section 2.1) in order to actually extract the concepts related to image processing. This extracting process is detailed more in-depth by the authors in [START_REF] Roynard | An Image Processing Library in Modern C++: Getting Simplicity and E ciency with Generic Programming[END_REF]. We performed the image processing concept extraction and made it available alongside the image processing library Pylene [START_REF] Carlinet | Pylena: a Modern C++ Library for Generic and E cient Image Processing[END_REF].

3 Image Views: Ease & Performance 3.1 Ranges and Views in the C++20 STL C++20 ranges [START_REF] Niebler | P1037R0: Deep Integration of the Ranges TS[END_REF] formalizes the concept of view, extending the array views implemented in array-manipulation libraries [START_REF] Andres | Runtime-Flexible Multi-Dimensional Arrays and Views for C++98 and C++0x[END_REF][START_REF] Veldhuizen | Blitz++: The Library that Thinks it is a Compiler[END_REF], and transferable to the Image concept. In the STL, there is a distinction between the container owning the data bu er, the iterators related to traversing this container, the range encapsulating the iterator pair allowing traversing the container and the view which mutates the way the base range traverse the data it is related to. All those abstraction levels need proper re ned design about data ownership and lifetime, depending on what it refers to. For instance, a range may not be cheap-to-copy as it may contain data to prolong lifetime of the underlying object, as for extending the lifetime of a temporary range in a pipe.

auto h = [](int x) 0 if x < 150 255 if x ≥ 150 auto u = auto v = transform( u , h ) → h u
In our design, all images have reference semantics and are cheap-to-copy. An image view, is a lightweight object that acts like an image and models the Image concept. For example, it can be a generator image object which generates a new value whenever ( ) is called, or an observer image that records the number of times each pixel is accessed to compare algorithms performance. In some pre C++-11 libraries (e.g. GIL [START_REF] Bourdev | Generic Image Library[END_REF] or Milena [START_REF] Levillain | Practical Genericity: Writing Image Processing Algorithms both Reusable and E cient[END_REF]), image views were also present (named morphers alongside the SCOOP pattern [START_REF] Burrus | A static C++ object-oriented programming (SCOOP) paradigm mixing bene ts of traditional OOP and generic programming[END_REF][START_REF] Géraud | Semantics-Driven Genericity: A Sequel to the Static C++ Object-Oriented Programming Paradigm (SCOOP 2)[END_REF]) but not compatible with modern C++ idioms (e.g. the rangebased for loop) and not as well-developed as in [START_REF] Niebler | P1037R0: Deep Integration of the Ranges TS[END_REF].

Among image views, we focus on image adaptors.

Let v = transform( 1 , 2 , • • • , Ĥ , ℎ)
where ğ are input images and ℎ a n-ary function. transform returns an image generated (adapting) from other image(s) as shown in g. 7. An adaptor does not "own" data but stores the transformation ℎ and references to the input images. The properties of the resulting view depend on ℎ. The projection h1: ( , , ) ↦ → that selects the green component of an RGB triplet gives a view that is writable and has the same domain as 1 . However, the projection h2: ( , ) ↦ → ( + )/2, applied on images 1 and 2 gives a read-only view that computes, pixel-wise, the average of 1 and 2 .

Following the same principle, a view can apply a restriction on an image domain. In g. 8, we show the adaptor clip(input, roi) that restricts the image to a non-regular roi and filter(input, predicate) that restricts the domain based on a predicate. All subsequent operations on those images will only a ect the selected pixels.

Views applied to image processing

Views feature many interesting properties that change the way we program. Let us consider the simple but real case image processing pipeline aiming at extracting objects from a background which is depicted in g. 9. Simply said, it computes the di erence between an image and a background image. The Gaussian blur and the morphological opening allow some robustness to noise. All the concepts and techniques shown here can be generalized to more complicated IP pipelines. Views are composable. One of the most important feature in a pipeline design (generally, in software engineering) is object composition. It enables composing simple blocks into complex ones. Those complex blocks can then be managed as if they were still simple blocks. In g. 9, we have 5 simple image operators Image → Image (grayscale conversion, subtract, thresholding, opening, Gaussian). As shown with MyComplexOperator, algorithm composition would consider these 5 simple operators as a single complex operator Image → Image that could then be used in another even more complex processing pipeline. Just like algorithms, image views are composable, e.g. a view of the view of an image is still an image. As shown with MyComplexImage1 & 2, we pipe the input images into other views before feeding them into algorithms .

Views improve usability. The code to compose images in g. 9 is as simple as what is shown in g. 10. People familiar with functional programming may notice similarities with these languages where transform (map) and lter are sequence operators. Views use the functional paradigm and are created by functions that take a function as argument: the operator or the predicate to apply pixel-wise; we do not iterate by hand on the image pixels. Also, it o ers the opportunity to raise the IP practitioner reasoning by one level. Indeed, in our example g. 10, the thresholding routine is written without any consideration of the dimension of the image (2D, 3D), nor the underlying data-type (8-bits, oat, RGB). Thanks to this new abstraction level, we are able to express more complex image processing routines very smoothly while remaining generic by default. For instance, restricting the input image to a speci c region or to a speci c color channel becomes as simple as writing the following code: auto imroi = thresholdf(view::clip(img_gray, roi)); // ROI auto imred = thresholdf(view::red(img_gray));//Red channel Views for lazy computing. Because the operation is recorded within the image view, this new image type allows fundamental image types to be mixed with algorithms. In g. 10, the creation of views does not involve any computation in itself but rather delays the computation until the expression ( ) is invoked (requesting the value for the point ). Because views can be composed, the evaluation can be delayed quite far. Image adaptors are template expressions [START_REF] Veldhuizen | Expression Templates[END_REF][START_REF] Veldhuizen | Blitz++: The Library that Thinks it is a Compiler[END_REF] as they record the expression used to generate the image as a template parameter. A view actually represents an expression tree (cf. g. 13). This approach is close to designing a Domain Speci c Language (DSL) [START_REF] Van Deursen | Domain-Speci c Languages: An Annotated Bibliography[END_REF], but it remains within the C++ languages, which is unlike another DSL for Image Processing named Halide [START_REF] Ragan-Kelley | Halide: A Language and Compiler for Optimizing Parallelism, Locality, and Recomputation in Image Processing Pipelines[END_REF] that have chosen to provide their toolchain infrastructure to solve performance issues related to heterogeneous computing.

Views for performance. With a classical design, each operation of the pipeline is implemented on "its own". Each operation requires memory to be allocated for the output image and also, each operation requires that the image is fully traversed. This design is simple, exible, composable, but is not memory e cient nor computation e cient. With the lazy evaluation approach, the image is traversed only once. That has three bene ts. First, there are no intermediate images, which is very memory e ective. Second, traversing the image maybe be faster thanks to a better memory cache usage. Indeed, a view acts as if we were writing an optimal operator that would combine all the operations. This approach is somewhat related to the kernel-fusing operations available in some HPC speci cations [START_REF]OpenVX[END_REF] but views-fusion is directly optimized by the C++ compiler [START_REF] Brown | Introducing Parallelism to the Ranges TS[END_REF]. Also, as we saw that the operations are represented by a functional Abstract Syntax Tree (AST), the compiler is able to perform optimization based on the deforestation algorithm [START_REF] Wadler | Deforestation: Transforming programs to eliminate trees[END_REF] that removes intermediate computation trees in order to achieve a treeless form. Third, and the most important, if there is a domain restricting operation downstream in a view pipeline, the operation upstream will only be performed on the relevant region of interest, instead of being performed on the whole image if there was no lazy-computing involved.

Reasoning at image level

The nal argument we bring in our discussion about views is the fact that the IP practitioner raises his reasoning by one level. Indeed, let us take a look at the alpha-blending algorithm as a support example for our argument. The default code for a classical, handmade (and error-prone C++) alphablending is presented in g. 11. This algorithm makes several non-relevant hypotheses about the image type. Indeed, it is not relevant to the nal application whether the image's color is 8-bits RGB or oat. Also, the practitioner may only need to process a speci c color channel, or a speci c region of the image. The image may also be 3D. To summarize, there are a lot of hypotheses that are not relevant to the application logic and yet weight on the resulting implementations which lead us to the need of genericity. The solution is to shift the abstract level by one layer and reason at image level, as shown in g. 12 which presents the code and the produced view expression tree. Rewriting the low level algorithm in terms of views is as simple as in g. [START_REF] Dehnert | Fundamentals of Generic Programming[END_REF]. Finally, we also show in g. 14 how simple it now becomes to restrict input images to a speci c region or speci c color channel directly by chaining views at image level thanks to being able to reason at image level. The code has become more readable, more expressive and more e cient by default.

Comparison with Data Flow oriented

frameworks A parallel can be drawn between image views and the data ow oriented programming [START_REF] Rey | MapReduce: Simpli ed Data Processing on Large Clusters[END_REF] style used in Data Science such as Apache Spark [START_REF] Zaharia | Resilient Distributed Datasets: A Fault-Tolerant Abstraction for In-Memory Cluster Computing[END_REF] or even TensorFlow [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems[END_REF]. Indeed, we nd similar properties in those data ow system, such as composition and lazy-computing. Data ow oriented programming is heavily inspired from functional programming collections (Scala, Haskell). Let us focus on the Apache Spark programming model for this comparison.

It is very similar to our view design in the fact that transformations (partitioning functions) can be compared to views (computed lazily and chainable) and actions (performed on Resilient Distributed Dataset (RDD) constructed from transformations) can be compared to our algorithms (perform the work and resolve the transformation). However, it differs from views at runtime. Views may only do computation on a speci c requested part of an image (as seen in section 3.2) whereas the Spark pipeline may do transformations on the whole dataset prior to performing a narrowing transformation. Indeed, Spark aims at distributing asynchronously computation on clusters, that is why it does not prevent inefcient and unnecessary computation due to the nature of the acyclic computation graph computed from RDDs. In contrast, views are static, and their compositions is static. There is no need for a framework for that. Also, Views enables in-place computation (through projectors) which is very memory e cient.

Finally, our design di ers in the sens that views are still image types (embedding an operation). The IP practitioner can focus on the behavior of his images and algorithms by reasoning at image level. On the other hand, data ow programmer focuses on the data and how to transform it in order to extract information. Design-wise, an RDD is a generalized super-type of data, more exible due to its dynamic nature, but it does not abstract away the underlying complexity incurred by the processed data.

Experimentation

To highlight the interest of Generic Programming and views in the context of performance-sensitive applications, we study the impact on a simple but real case image processing pipeline aiming at extracting objects from a background as depicted on g. 9. Simply said, it computes the di erence Table 2. Benchmarks of the pipeline g. 9 on a dataset (12 images) of 10MPix images with OpenCV as a baseline.

between an image and a registered image. The Gaussian blur and the morphological opening allow some robustness to noise. The pipeline is implemented with (1) OpenCV, (2) our library (Pylene) where each step is a computing operator, (3) and our library where the purple blocks are views. This pipeline actually produces interesting results, as shown in g. 15. In table 2, we benchmark the computation time and the memory usage (using valgring/massif ) of these implementations (all single-threaded) with an opening of disc of radius 32 on 10 MPix RGB images (the minimum of many runs is kept). The results should not be misunderstood. They do not say that OpenCV is faster or slower but shows that implementations all have the same order of processing time, despite the fact that the algorithms used in our implementation are not the same as those used in OpenCV for blur and dilation/erosion, so that the comparison makes sense. It allows us to validate experimentally the advantages of views in pipelines. First, we have to be cautious about the real bene t in terms of processing time. Here, most of the time is spent in algorithms that are not eligible for view transformation. Thus, depending on the operations of the pipeline, views may not improve processing time. Nevertheless, using views does not degrade performance neither (only 1% in this experiment). It seems to show that using views does not introduce performance penalties and may even be bene cial in lightweight pipelines as the one in section 3. On the memory side, views reduce drastically the memory usage which is bene cial when developing applications which are memory constrained. On the binary size side, while we can fear "code bloat" due to many template instanciations, it appears that our binary size is very much in the same order of magnitude as for one doing the same work with OpenCV. We can deduce that the compiler is able to prune generated code and keep only needed and optimized machine code in the nal binary. From the developer standpoint, it requires only few changes in the code as shown in g. 10 -the implementation of the algorithms remains the same -which is a real advantage for software maintenance.

To conclude, the main advantage of views lies in the usability gain for the user without degrading performance. They can even improve memory usage e ciency, which is useful for embedded systems. To showcase this point, we o er additional material 1 demonstrating how views can, for instance, be used for a connected-component labeling algorithm. On a toy example, (a simpler pipeline) we show that views can improve performance (the source is in the supplementary material).

Limitations & Conclusion

This paper demonstrates the usage of image views, a modern C++ feature applied to image processing, that reconciles genericity, e ciency and ease of use. This contribution allows the IP practitioner to reason at image level, which is novel in the C++ IP library ecosystem, and facilitates the transition between prototypes and production-ready applications. In addition to improving the usability, it leads to performance gain and decreases the memory usage. These ideas have been implemented in our C++20 library [START_REF] Carlinet | Pylena: a Modern C++ Library for Generic and E cient Image Processing[END_REF] and used for concrete image processing applications (medical imaging and document analysis).

Despite these advantages, this design comes with some limitations. The intensive use of the C++ templates generates type combinatorial explosion and code bloat that leads to large compilation times. Also, templates belong to the static world and are poorly interoperable with the dynamic world. Nevertheless, dealing with the dynamic world (runtime) is mandatory when it comes down to exposing a static library to a dynamic language like Python, in order to provide the interactivity that C++ lacks. Another drawback lies in relying on the compiler to do a lot of work, especially the optimizer. Therefore, when the compiler fails in optimizing our abstractions away, we pay the price in performance. Our future work will be dedicated to x these shortcomings.

Figure 1 .

 1 Figure 1. The watershed segmentation algorithm runs on a 2D-regular grayscale image (left), on a vertex-valued graph (middle) and on a 3D mesh (right).

Figure 4 .Figure 5 .

 45 Figure 4. Concept-checked sum algorithm over a collection.

Figure 6 .

 6 Figure 6. Image concept and constrained dilation algorithm.

( p o i n t e r ) ≡Figure 7 .Figure 8 .

 )78 Figure 7. An image view performing a thresholding.

Figure 9 .Figure 10 .

 910 Figure 9. Pipeline for foreground extraction using algorithms and views. float kThreshold = 150; float kVSigma = 10; float kHSigma = 10; int kOpeningRadius = 32; auto img_gray = view::transform(img_color, to_gray); auto bg_gray = view::transform(bg_color, to_gray); auto bg_blurred = gaussian2d(bg_gray, kHSigma, kVSigma); auto tmp_gray = img_gray -bg_blurred; auto thresholdf = [](auto x) { return x < kThreshold; }; auto tmp_bin = view::transform(tmp_gray, thresholdf); auto ero = erosion(tmp_bin, disc(kOpeningRadius)); dilation(ero, disc(kOpeningRadius), output);Figure 10. Pipeline implementation with views .

Figure 11 .Figure 12 .Figure 13 .

 111213 Figure 11. Alpha-blending with classical C/C++ code.

Figure 14 .Figure 15 .

 1415 Figure 14. Chaining views to feed alpha-blending.

  , let us consider the morphological dilation of an image : → (de ned on a domain with values in ) by

	template <class I>
	concept Image = requires {
	point_t<I>; // Type of point (P)
	value_t<I>; // Type of value (V)
	} && requires (I f, point_t<I> p, value_t<I> v) {
	{ v = f(p) }; //
	{ f(p) = v }; // optional, for output
	{ f.domain() } -> Range; // (actually Range of P)
	};
	template <Image I, StructuringElement SE>
	void dilation(I input, I output, SE se)
	requires MaxMonoid<value_t<I>>
	{ ... }

December 06-07, 2022, Auckland, New Zealand. ACM, New York, NY, USA, 8 pages. h ps://doi.org/10.1145/3564719.3568692


Baptiste Esteban

Guillaume Tochon

Edwin Carlinet

Didier Verna

Estimation of the noise level function for color images using mathematical morphology and non-parametric statistics

Noise level information is crucial for many image processing tasks, such as image denoising. To estimate it, it is necessary to find homegeneous areas within the image which contain only noise. Rank-based methods have proven to be efficient to achieve such a task. In the past, we proposed a method to estimate the noise level function (NLF) of grayscale images using the tree of shapes (ToS). This method, relying on the connected components extracted from the ToS computed on the noisy image, had the advantage of being adapted to the image content, which is not the case when using square blocks, but is still restricted to grayscale images. In this paper, we extend our ToS-based method to color images. Unlike grayscale images, the pixel values in multivariate images do not have a natural order relationship, which is a well-known issue when working with mathematical morphology and rank statistics. We propose to use the multivariate ToS to retrieve homogeneous regions. We derive an order relationship for the multivariate pixel values thanks to a complete lattice learning strategy and use it to compute the rank statistics. The obtained multivariate NLF is composed of one NLF per channel. The performance of the proposed method is compared with the one obtained using square blocks, and validates the soundness of the multivariate ToS structure for this task.

I. INTRODUCTION

During acquisition, natural images may be corrupted by some noise due to physical phenomena. Knowledge of the parameters of this noise related to the image pixels intensity, called noise level, is important for image processing applications such as denoising [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], segmentation [START_REF] Droske | Multiscale joint segmentation and registration of image morphology[END_REF], super-resolution [START_REF] Freeman | Example-based superresolution[END_REF] or image compression [START_REF] Walker | Combined image compressor and denoiser based on treeadapted wavelet shrinkage[END_REF] and can greatly improve the performance of these algorithms. Furthermore, knowledge of the statistical nature of the noise is also important as it changes the way these applications are processed.

In the past, several methods to estimate the noise level have been proposed. Some of them separate the noise from the image content, either in the case of a signal-independent noise [START_REF] Liu | Single-image noise level estimation for blind denoising[END_REF] or a signal-dependent noise such as the Poisson-Gaussian noise [START_REF] Mäkitalo | Noise parameter mismatch in variance stabilization, with an application to poisson-gaussian noise estimation[END_REF], using a variance stabilization method. With the increasing usage of deep learning based methods for image processing algorithms, some noise level estimations relying on convolutional neural networks were recently proposed, usually in the context of noise removal, for additive Gaussian noise either independent and identically distributed [START_REF] Zhang | Ffdnet: Toward a fast and flexible solution for cnn-based image denoising[END_REF], [START_REF] Xu | A fast yet reliable noise level estimation algorithm using shallow cnn-based noise separator and bp network[END_REF] or not [START_REF] Yue | Variational denoising network: Toward blind noise modeling and removal[END_REF], [START_REF] Tan | Pixelwise estimation of signal-dependent image noise using deep residual learning[END_REF]. Even though these methods are very effective, they still require knowledge of the statistical nature of the noise beforehand. Beaurepaire et al. [START_REF] Beaurepaire | Identification of the nature of noise and estimation of its statistical parameters by analysis of local histograms[END_REF] propose to first identify the nature of the noise and then estimate its level by analizing statistics on patches of image, but it did not allow to represent a mixed noise. Finally, Sutour et al. [START_REF] Sutour | Estimation of the noise level function based on a nonparametric detection of homogeneous image regions[END_REF] propose a method to estimate a noise level function (NLF) [START_REF] Liu | Noise estimation from a single image[END_REF], represented as a quadratic function, from statistics computed from homogeneous square patches detected by a non-parametric statistical test, which is extended by Esteban et al. [START_REF] Esteban | Estimating the noise level function with the tree of shapes and non-parametric statistics[END_REF] to homogeneous regions to fit the image content. However, this method is still limited to univariate images.

In this paper, we propose to extend our previous method to color images, solving the issues raised by the multivariate pixel values such as the lack of natural order relationship using a rank map resulting from a learned complete lattice [START_REF] Lézoray | Complete lattice learning for multivariate mathematical morphology[END_REF] and the usage of the multivariate tree of shapes [START_REF] Carlinet | Mtos: A tree of shapes for multivariate images[END_REF] for the image simplification. Furthermore, we improve the detection of homogeneous regions, using the adjacency relationship as in [START_REF] Sutour | Estimation of the noise level function based on a nonparametric detection of homogeneous image regions[END_REF]. Then, we define a noise level function for color images constrained to the case where the noise is not correlated between the channels. Finally, we evaluate the performance of our proposed method on a benchmark of 150 natural color images.

Thus, this paper is structured as follow: first, we recall the estimation of the NLF for univariate images in Section II. Then, its extension to color images is described in Section III. We measure the performance of this extension by comparing it with different methods in Section IV. Finally, we conclude in Section V.

II. ESTIMATION OF THE NLF FOR UNIVARIATE IMAGES

This section recalls the main steps of [START_REF] Esteban | Estimating the noise level function with the tree of shapes and non-parametric statistics[END_REF] for the regionbased NLF estimation using the tree of shapes. The noise level of an image represents the variance of the noise related to the pixel intensities of this image. It may be represented as a function, the noise level function (NLF) [START_REF] Liu | Noise estimation from a single image[END_REF]. This function is modeled as a quadratic function by Sutour et al. [START_REF] Sutour | Estimation of the noise level function based on a nonparametric detection of homogeneous image regions[END_REF], where each coefficient corresponds to the variance of one type of statistical noise.

Let θ = (a, b, c) ∈ (R + ) 3 . The noise level function of an image is defined as

σ 2 = NLF θ (x) = ax 2 + bx + c (1)

A. The tree of shapes

The tree of shapes (ToS) [START_REF] Monasse | Fast computation of a contrast-invariant image representation[END_REF], denoted by T , is a hierarchical representation encoding the inclusion relationship of the connected components of an image. Let f : Ω → R be an image defined on a domain Ω ⊂ Z 2 . The upper level set is defined as [f ≥ λ] = {x ∈ Ω|f (x) ≥ λ} and the lower level set as [f < λ] = {x ∈ Ω|f (x) < λ}. A shape C ⊂ Ω is a connected component belonging to one of these sets with its holes filled. Fig. 1b illustrates the ToS of the image Fig. 1a.

Furthermore, the contours of the shapes, denoted by ∂C, represent the level lines in the image. Thus, the ToS also encodes the inclusion relationship of the level lines of an image. Finally, the ToS is computed with a linear complexity [START_REF] Carlinet | The tree of shapes turned into a max-tree: a simple and efficient linear algorithm[END_REF].

B. Simplification of the ToS

The meaningfulness of a shape can be its importance in the description of the image. A meaningfulness measure is the average gradient magnitude along the contour of a shape, which represents a level line of the image. Thus, removing meaningless shapes from the ToS is a simplification process.

A mean to simplify the ToS rely on the piecewise-constant Mumford-Shah functional [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF]. It is defined as below:

E(R) = Ri∈R x∈Ri ∥ f (R i ) -f (x)∥ 2 2 + λ|∂R| (2)
with R = ⊔R i for R i a connected component in the domain of f , f (R i ) is the mean of the values inside a region R i , ∂R = ∪∂R i is the union of the boundaries of each region R i and λ is a regularization parameter.

When the minimization of this functional is subordinated to the ToS, the process of simplification consists in selecting a subset of level lines. Thus, the regions boundaries are then fixed to the level lines ∂C i . Formally, the minimization process is given by: min

Ts E(T s ) (3)
where T s is the simplified ToS. This simplification process is illustrated Fig. 1c, where the blue dashed lines are related to the parenthood relationship before the simplification and the plain blue line to this relationship after. The resulting reconstructed image Fig. 1d illustrates the impact of the removal of the node C in T . An algorithm has been proposed by Xu et al. [START_REF] Xu | Salient level lines selection using the mumford-shah functional[END_REF] to simplify the ToS, sorting the level lines in a meaningful order and iteratively removing the shapes to decrease the energy of the functional.

C. Detection of homogeneous regions

Once the ToS has been simplified, the connected components whose size is above a given threshold are extracted and their homogeneity is tested. A region is homogeneous if it contains only noise, meaning that pixel values are not correlated between them. This correlation is measured using the Kendall τ rank correlation coefficient [START_REF] Kendall | A new measure of rank correlation[END_REF].

Let (x, y) ∈ (R n ) 2 be two sequences of observations of size n of two random variables X and Y and x i ∈ x is the i th value of x and y i ∈ y be the i th value of y. The Kendall τ rank correlation coefficient is defined on the interval [-1, 1] by:

τ (x, y) = 1 n -1 1≤i,j≤n sign(x i -x j ) sign(y i -y j ) (4)
for x i ̸ = x j and y i ̸ = y j . Furthermore, the definition of the Kendall τ has been extended to take into account tied pairs (x i = x j or y i = y j ) [START_REF]The treatment of ties in ranking problems[END_REF]. The detection of homogeneous regions relies on a statistical test associated with the Kendall τ . Under the null hypothesis H 0 , stating the absence of correlation between the two sequences x and y (τ (x, y) = 0), the z-score associated to the Kendall τ follows a standard normal distribution N (0, 1). In contrast with H 0 , the alternative hypothesis H 1 is declared when τ (x, y) ̸ = 0, indicating the fluctuation of the values cannot be explained by the noise only.

To perform the detection, a region is divided randomly into two sequences x and y of n pixels several times. The p-value p = P(τ (X, Y ) > τ (x, y) | H 0 ) is computed for all divisions and if the second smallest p-value is superior to a predefined level of detection α, the region is considered homogeneous. In practice, each region is divided randomly 10 times.

D. Estimation of the NLF

The NLF estimation is a 3 steps process: the segmentation of the image, the statistical test on the connected components extracted from these segmentations and the estimation of the NLF using statistics computed from the homogeneous regions.

The segmentation process is as follow: first, the image is blurred with a gaussian kernel, reducing the amount of level lines in the image. Then, a ToS is built on this blurred image. This ToS is simplified using the Mumford-Shah functional with various values of the regularization parameter λ. This allows to obtain image segmentations of various simplification scales, with regions of different size. The number of removed shapes in the tree increases with the value of the regularization parameter λ, resulting in a more simplified reconstructed image.

Then, the homogeneity of the regions extracted from the simplified ToS is tested using the Kendall τ , but using the pixel values of the noisy image to conduct the rank-based test. The first and second order statistics, denoted respectively by μ and σ2 of those homogeneous regions are computed, leading to information about the noise. From Eq. 1, the NLF is the variance of the noise related to pixel intensity. Thus, the NLF is estimated by minimizing the residual error between the empirical variance of a region σ2

i and the variance predicted by the relation NLF θ ( μi ) for all homogeneous region:

θ = (a, b, c) = argmin θ∈(R + ) 3 i ∥NLF θ ( μi ) -σ2 i ∥ 1 ( 5
)
The minimization is performed using the preconditionned primal-dual algorithm of Chambolle-Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], relying on a L 1 minimization since it is more robust to outliers than a L 2 minimization.

Comparative results between the block-based estimation from Sutour et al. and the estimation described in this section were presented in [START_REF] Esteban | Estimating the noise level function with the tree of shapes and non-parametric statistics[END_REF] and demonstrated the pertinence in using regions adapted to the image content instead of square blocks.

III. EXTENSION OF THE ESTIMATION TO COLOR IMAGES

Color images suffer from the multivariate nature of their values. Rank-based statistics and mathematical morphology are two fields strongly relying on the rank of the considered values. However, multivariate values do not have a natural ordering, which is a well known and studied issue in mathematical morphology [START_REF] Aptoula | A comparative study on multivariate mathematical morphology[END_REF]. In order to make it possible to estimate the NLF for color images, we use tools from multivariate mathematical morphology: first, the ToS has been extended to multivariate images, resulting in the multivariate ToS [START_REF] Carlinet | Mtos: A tree of shapes for multivariate images[END_REF], its construction depending on the inclusion relationship of the shapes obtained from the channelwise ToSs. Then, we propose to construct a total ordering on the pixel color values using a complete lattice learning approach [START_REF] Lézoray | Complete lattice learning for multivariate mathematical morphology[END_REF], which allows us to extend the computation of the Kendall τ rank correlation coefficient.

A. Rank transform based on complete lattice learning

Let f : Ω → I ⊂ R n be a multivariate image defined on a domain Ω and I its pixel values set containing m vectors of n element in R. A rank transform R : I → [1, m] is a mapping from an element of I to its rank according to the order relationship endowing I. This operation requires I to Lezoray [START_REF] Lézoray | Complete lattice learning for multivariate mathematical morphology[END_REF] proposes to solve this issue by using an unsupervised h-ordering. A h-ordering, denoted by ≤ h , is the application of a bijective projection h : I → L such that, for any two elements u, v ∈ I, u ≤ h v ⇔ h(u) ≤ h(v). Thus, I does not need to be a complete lattice anymore, since the ordering relationship depends on the elements of L.

The construction of the projection h is a three steps process. First, a dictionary D of p ≪ m elements is extracted out of the whole amount of data using a vector quantization algorithm (such as the k-means clustering algorithm for instance). This quantization is necessary due to the complexity of the manifold learning algorithm, which is quadratic with the number of input samples.

The new representation is constructed through dimensionality reduction thanks to the Laplacian Eigenmaps [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF]. To this aim, the pairwise similarity between each element of the dictionary is computed using a Gaussian kernel, resulting in a similarity matrix S D . The similarity matrix can be seen as the adjacency matrix of a graph, where each node is an element of D. The diagonal degree matrix, denoted by D D , is the sum of the weights of the edges connected to each node such that (D D ) ii = j (S D ) ij . Using these matrices, the Laplacian matrix is obtained computing

L D = D D -S D .
The new representation is finally obtained by computing the eigenvectors of the normalized Laplacian matrix defined as follows:

L D = D -1 2 D L D D -1 2 D (6)
The obtained representation, denoted by h D , is defined by

h D (c i ) = (ϕ 1 D (c i ), ..., ϕ p D (c i )) where c i is the i th element Input image T1 T2 T3 GoS G 0 1 2 2 1 3

Depth computation

Depth map

Hole-filled max-tree Fig. 3: Workflow of the MToS construction of D and ϕ j D is the j th eigenvector of L D , with ϕ j D (c i ) its value at index i. However, it is only valid for the dictionary values resulting from the vector quantization. Consequently, it needs to be extended to all image values. This is done using the Nyström extension [START_REF] Talwalkar | Large scale svd and manifold learning[END_REF]. For Φ D the matrix of the eigenvectors, Λ D the diagonal matrix of eigenvalues obtained by the eigen-decomposition of the normalized Laplacian matrix and D I is the diagonal degree matrix computed from I, the Nyström extension is defined by:

Φ I = D -1 2 I S T D D -1 2 D Φ D Λ -1 D (7)
The projection h is defined as h(v) = (ϕ 1 I (v), ..., ϕ p I (v)) where v ∈ I and ϕ i I is the i th extrapolated eigenvector of the normalized Laplacian matrix for all the values in I. Finally, the complete lattice (I, ≤) is obtained by using the lexicographic ordering on this new representation and the rank map R f resulting from the application of the rank transform on the image f is computed. The results of each step of the complete lattice learning to obtain the rank map are illustrated in Fig. 2. Instead of computing the Kendall τ directly on the image value, this rank map is used.

B. Multivariate Tree of Shapes

Instead of using the total ordering provided by the rank map, the ToS is extended to multivariate images by using the inclusion relationship between the shapes provided by the channelwise ToSs. This new representation, proposed by Carlinet et al. [START_REF] Carlinet | Mtos: A tree of shapes for multivariate images[END_REF], is the multivariate tree of shapes (MToS). Its computation, described below, is composed of two major steps: the computation of a graph of inclusion and the extraction of a new tree from this graph.

First, the ToS T i of each channel f i is computed, with C i its associated set of shapes. Then, these trees are merged into a graph of shapes G, which is a directed acyclic graph. Its edges represent the inclusion relationship between the connected components of a given channel as for the ToS, but also the inclusion of the shapes between the channelwise ToSs. The root of the trees, representing the whole image, constitutes the root of this graph. Thus, for C = i C i , the set (C, ⊆) is the cover of G.

The second step consists in extracting a tree from G. First, the depth attribute is computed. As G is a directed acyclic graph, this attribute corresponds to the longest path from a

(a) (b) (c) (d) (e)
Fig. 4: Region division schemes into two sequences x and y node to the root of G. Then, a depth map is built, associating to each point of the original image the maximum depth value of the nodes associated to the current point. Finally, a maxtree is built on the depth map, ensuring that the components of the upper level set of the map are valid shapes by applying a hole-filling operator. The resulting tree is MToS. The MToS computation is summarized in Fig. 3.

C. Estimation of the color noise level function

The NLF defined by ( 1) is a valid representation for univariate images, but does not allow to represent the noise level for color images. Let CNLF(x) = (NLF θr (x r ), NLF θg (x g ), NLF θ b (x b )) be the color noise level function (CNLF) where x ∈ R 3 is a multivariate pixel value and subscript r, g, and b correspond respectively to the red, green and blue channels of the color image. This representation assumes that the noise is not correlated between the channels.

The tools described previously are used to extend the estimation of the NLF to color images. First, the MToS is used instead of the ToS. It provides the simplifications to extract the connected components. Then, each connected component is divided into two sequences using the adjacency relationship scheme proposed in Sutour et al. [START_REF] Sutour | Estimation of the noise level function based on a nonparametric detection of homogeneous image regions[END_REF]. Thus, the division is processed 4 times, using the adjacent relationship as illustrated by Fig. 4a to 4d for the horizontal, vertical and the two diagonal directions, and the Kendall τ rank correlation coefficient is computed using the rank map obtained from a rank transform on the original image.

Finally, the empirical mean and variance are computed on each channel for the detected homogeneous regions and the minimization problem expressed by equation ( 5) is solved channelwise using the Chambolle-Pock algorithm.

IV. RESULTS

A. Experimental set-up

Following what is done in [START_REF] Sutour | Estimation of the noise level function based on a nonparametric detection of homogeneous image regions[END_REF] and [START_REF] Esteban | Estimating the noise level function with the tree of shapes and non-parametric statistics[END_REF], we use the mean relative error (MRE) to compare two CNLF channelwise in the following experiments. Let NLF θ be the NLF corresponding to a known set of parameters θ = (a, b, c) and NLF θ be the estimated NLF. The MRE is defined as:

MRE( θ) = 1 |I| x∈I |NLF θ (x) -NLF θ (x)| NLF θ (x) ( 8
)
where I is the set of pixel values of the image. The lower the MRE, the better the estimation. The experiments are performed on all 150 images from the Laurent Condat image database 1 , where the acquisition noise in the high definition images can be neglected, making this dataset ideal to experiment on noise estimation tasks.

In the following experiments, the images from the benchmark data set are first polluted by a mixed noise with a given CNLF where the coefficients θ i = (a i , b i , c i ) are picked from a standard normal law N (0.01, 0.003). Following what is done for grayscale images in [START_REF] Sutour | Estimation of the noise level function based on a nonparametric detection of homogeneous image regions[END_REF], [START_REF] Esteban | Estimating the noise level function with the tree of shapes and non-parametric statistics[END_REF], the minimum size of the tested regions is set to 250 and the detection level α to 0.4. The MToS is simplified with a Mumford-Shah regularization parameter λ empirically set to 200, 500 and 700. To evaluate in practice the robustness of the proposed method, this experiment (adding a mixed noise, estimating the NLF, measuring the MRE) is repeated 20 times per image in the benchmark, for 20 different drawings of the noise parameters (a, b, c), and the average MRE for each channel of the image is computed. In order to present concisely and qualitatively the results obtained on the whole benchmark, we use 3 scatterplots (on per channel) where x-axis and y-axis support the average MRE for the two methods being compared. Thus, a point above the plotted blue line corresponds to better results for the method supported by the x-axis in the scatterplot.

B. Evaluation of the division scheme

In [START_REF] Sutour | Estimation of the noise level function based on a nonparametric detection of homogeneous image regions[END_REF], Sutour et al. propose to use the adjacency schemes to perform the division of blocks into sequences x and y. Because the regions extracted from the ToS do not have regular shapes allowing these division schemes, we propose in [START_REF] Esteban | Estimating the noise level function with the tree of shapes and non-parametric statistics[END_REF] to perform a random splitting instead, illustrated by Fig. 4e.

In the current work, we compare the random splitting of [START_REF] Esteban | Estimating the noise level function with the tree of shapes and non-parametric statistics[END_REF] with the adjacency-based division of [START_REF] Sutour | Estimation of the noise level function based on a nonparametric detection of homogeneous image regions[END_REF], adapted to the irregular morphology of the shapes extracted from the MToS. Even though the random division gave better results than the adjacency-based division of blocks for grayscale images, the adjacency relationship improves the accuracy of the estimation, according to the average channelwise MREs. Thus, for the red, green and blue CNLF components, the average MRE of the adjacency-based division scheme is better for respectively 146, 149 and 147 images.

C. Comparison with a block-based estimation

In order to be concordant with the experiments performed in [START_REF] Esteban | Estimating the noise level function with the tree of shapes and non-parametric statistics[END_REF], we evaluate the accuracy of the proposed CNLF estimation with an extension of the block-based estimation proposed in [START_REF] Sutour | Estimation of the noise level function based on a nonparametric detection of homogeneous image regions[END_REF]. The difference between the two compared methods is illustrated in Fig. 6. The block-based estimation divides the image into several square blocks. Each block is then split into two sequences of pixels using an adjacencybased scheme. The statistical test is then performed on these sequences, but instead of using the rank map resulting from a rank transform, the Kendall τ rank correlation coefficient is computed on the pixel intensity values of each channel, since it turned out to be more efficient for the block-based estimation. A block is declared homogeneous if all the three channelwise statistical tests are concordant and in favor of homogeneity. In order to attain an overall detection level of α, this detection level is lowered to α/3 for the three channelwise tests. Futhermore, we extended the automatic block size decrease by imposing the same criteria as in [START_REF] Sutour | Estimation of the noise level function based on a nonparametric detection of homogeneous image regions[END_REF] for every channel.

Quantitative results between the region-based estimation and the block-based estimation are shown in Fig. 7. For the red, green and blue CNLF components, the average MRE is better for respectively 147, 137 and 139 images. These results confirm that using the MToS to adapt the detection of homogeneous areas to the various content of the 150 benchmark images improves the performance of the estimation compared to the usage of blocks.

Furthermore, the difference between the estimations for one image is illustrated in Fig. 8. The homogeneous regions shown in Fig. 8b are extracted from all simplifications. In the plots of the resulting CNLF in Fig. 8c and8e, there is less statistical samples (19 samples) used to estimate the channelwise NLF for the region-based estimation than for the block-based estimation (46 samples) but they are more precise and the estimated channelwise NLFs, in plain lines in the plot, are closer to the reference NLFs, in dashed lines, used to add noise to the image. The MREs obtained with the proposed estimation are for each channel respectively 0.036, 0.039 and 0.048, which are better than the blockbased estimation's MREs, which are 0.118, 0.067 and 0.115. This confirms that using content-adapted homogeneous areas is more efficient than using blocks, even though the number of statistical samples is fewer. The obtained regions are larger than the fixed size blocks, resulting in more precise statistical means and variances used to estimate the CNLF.

The computational time for the block-based estimation ranges between 1.3 and 102.4 seconds, depending whether the block size was automatically decreased or not. For our proposed method, the computational time in around 66.4 seconds in total, decomposed as 10.4 seconds for the MToS computation and simplifications, 41.8 seconds for the complete lattice learning resulting in the rank map, 14.2 seconds for the homogeneous regions detection and 0.007 seconds for the minimization of Eq. 5.

V. CONCLUSION

In this article, we proposed to extend our previous method to estimate the NLF with the ToS and rank statistics to color images. Issues were raised by the multivariate nature of color However, this method extracts a partition of the tree which is optimal in the sense of the Mumford-Shah functional, but this partition depends on the value of a regularization parameter which greatly impacts the MToS simplification. Thus, only a subset of the regions contained in the MToS are effectively tested. Consequently, in our future work, we will extend the search for homogeneous areas to all regions contained in the tree structure. Furthermore, we will consider another kind of hierarchical structure being hierarchy of partition, whose construction relies on dissimilarity measures between neighborhing regions (and not on the absolute values of their pixels), and which lend themselves better to multivariate data. As the noise is often correlated between channels in real image processing applications, we plan to extend our CNLF model to correlated noise between channel and estimate it using the detected homogeneous regions, these ones handling the values of the three channels.

Fig. 1 :

 1 Fig. 1: Illustration of the simplification of the ToS

Fig. 2 :

 2 Fig. 2: Illustration of the complete lattice learning to obtain the rank map R f

Fig. 5 :

 5 Fig. 5: Comparison of the MRE of the shape-based estimation using a random division and an adjacency-based division For a matter of style, one point has been removed from the plots, respectively at (5.7, 2.2), (7.8, 10.4) and (6.1, 8.0).

Fig. 6 :

 6 Fig. 6: Summary of the estimation pipeline

Fig. 7 :

 7 Fig. 7: Comparison of the MRE between the block-based channelwise estimation and the region-based estimation. For a matter of style, one point has been removed from the plots, respectively at (3.9, 0.9), (5.0, 1.0) and (2.1, 1.0).

Fig. 8 :

 8 Fig. 8: Illustration of the region-based estimation and the block-based estimation