Markus Gotz
email: m.goetz@fz-juelich.de

Member, IEEE Gabriele Cavallaro
email: g.cavallaro@fz-juelich.de

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

Matthias Book
email: book@hi.is

Member, IEEE Morris Riedel
email: m.riedel@fz-juelich.de

Markus G Ötz

Member, IEEE Thierry G Éraud

M Götz

Parallel Computation of Component Trees on Distributed Memory Machines

Keywords: Component-Trees, Threshold Decomposition, Max-Tree, Connected Component Labeling, High-Performance Computing, Hybrid Application, MPI, Multithreading. ✦

HAL is

INTRODUCTION

S INCE the 1960s, mathematical morphology [1], [2] has become increasingly popular in the image processing community mainly due to its proven utility and rigorous mathematical description. The mathematical morphology framework provides a set of powerful operators for analyzing the spatial domain of images at the region-leveli.e., connected components-based on tree representations, called thresholds decompositions [START_REF] Breen | Attribute Openings, Thinnings, and Granulometries[END_REF], [START_REF] Salembier | Antiextensive Connected Operators for Image and Sequence Processing[END_REF]. These are based on tree representations of images which can be divided into two main groups [START_REF] Najman | A Graph-based Mathematical Morphology Reader[END_REF]: hierarchies of segmentation-i.e., hierarchy of image partitions such as minimum spanning trees [START_REF]On the Shortest Spanning Subtree of a Graph and the Traveling Salesman Problem[END_REF], alpha-trees [START_REF] Ouzounis | The Alpha-Tree Algorithm[END_REF], binary partition trees [START_REF] Salembier | Binary Partition Tree as an Efficient Representation for Image Processing, Segmentation, and Information Retrieval[END_REF]-and threshold decompositions-i.e., hierarchy of regions such as component trees [START_REF] Salembier | Antiextensive Connected Operators for Image and Sequence Processing[END_REF], [START_REF] Jones | Component Trees for Image Filtering and Segmentation[END_REF], tree of shapes (ToS) [START_REF] Caselles | Topographic Maps and Local Contrast Changes in Natural Images[END_REF] and multivariate tree of shapes [START_REF] Carlinet | MToS: A Tree of Shapes for Multivariate Images[END_REF]). Generally, tree structures are often considered richer in descriptive ability since they can be exploited for breaking down images into their fundamental elements which are easier to interpret with regards to the pixels. Component trees [START_REF] Salembier | Antiextensive Connected Operators for Image and Sequence Processing[END_REF], [START_REF] Jones | Component Trees for Image Filtering and Segmentation[END_REF], are thresholds decompositions that represent connected components [START_REF] Salembier | Flat Zones Filtering, Connected Operators, and Filters by Reconstruction[END_REF] at every threshold level of an image in a hierarchical fashion, through parent relationships between nodes. The connected components organized in such trees can be filtered with different strategies [START_REF] Breen | Attribute Openings, Thinnings, and Granulometries[END_REF], [START_REF] Salembier | Antiextensive Connected Operators for Image and Sequence Processing[END_REF] and can model various types of connectivity [START_REF] Ouzounis | Mask-Based Second-Generation Connectivity and Attribute Filters[END_REF].

Component trees (i.e., max-and min-tree) have been popularized by connected operators, such as attribute filters [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF], [START_REF] Breen | Attribute Openings, Thinnings, and Granulometries[END_REF], which have been extensively used for the modeling of spatial information in images from remote sensing [START_REF] Dalla Mura | Morphological Attribute Filters for the Analysis of Very High Resolution Remote Sensing Images[END_REF], [START_REF] Benediktsson | Hierarchical Analysis of Remote Sensing Data: Morphological Attribute Profiles and Binary Partition Trees[END_REF], astronomy [START_REF] Berger | Effective Component Tree Computation with Application to Pattern Recognition in Astronomical Imaging[END_REF], [START_REF] Teeninga | Improved Detection of Faint Extended Astronomical Objects Through Statistical Attribute Filtering[END_REF] and medical scanning [START_REF] Purnama | Non-Compactness Attribute Filtering to Extract Retinal Blood Vessels in Fundus Images[END_REF], [START_REF] Kiwanuka | Automatic Attribute Threshold Selection for Morphological Connected Attribute Filters[END_REF]. Attribute filters are edge-preserving and flexible operators due to the preservation of contours in the processed objects and rely on multiple spatial measures or attributes. The possibility to perform a multi-attribute analysis, like attribute filters built by employing different attributes, enriches the extraction of spatial arrangement and improves the discrimination between different structures. In the presence of scenes with high complexity and heterogeneity, e.g., densely populated urban area, a complete modeling of the spatial information can be achieved through a multi-level analysis. It implies the decomposition of the original gray-level image obtained by applying a sequence of attribute filters according to a set of filter thresholds [START_REF] Cavallaro | Automatic Attribute Profiles[END_REF]. The result of this operation are the so-called attribute profiles [START_REF] Dalla Mura | Morphological Attribute Filters for the Analysis of Very High Resolution Remote Sensing Images[END_REF]. They have been exploited mainly in remote sensing, e.g., classification [START_REF] Song | Remotely Sensed Image Classification Using Sparse Representations of Morphological Attribute Profiles[END_REF], [START_REF] Falco | Spectral and Spatial Classification of Hyperspectral Images Based on ICA and Reduced Morphological Attribute Profiles[END_REF], data fusion [START_REF] Pedergnana | Classification of Remote Sensing Optical and LiDAR Data Using Extended Attribute Profiles[END_REF] and change detection [START_REF] Falco | Change Detection in VHR Images Based on Morphological Attribute Profiles[END_REF], as well as in medical imaging processing for tomographic image segmentation [START_REF] Ouzounis | Differential Area Profiles: Decomposition Properties and Efficient Computation[END_REF]. Recently, attribute filters are utilized within a novel deep learning framework for the large-scale, unsupervised detection of objects in remote sensing image [START_REF] Aldeborgh | Unsupervised object detection on remote sensing imagery using hierarchical image representations and deep learning[END_REF]. A set of attributes is automatically identified in order to extract a representative, high quality training data set.

Nowadays, image processing applications rely on very high resolution data due to the continuing technological improvements of the sensor instruments. For example earth observation platforms have led to the increasing volume, acquisition speed and variety of sensed images, e.g., the (c) The corresponding max-tree. Fig. 1: Example of max-tree representation based on exemplary image and its components C i c , with the subscript c being the gray-level and the superscript i the canonical point uniquely identifying the component.

World-View-3 satellite sensor (spatial resolution of 0.31m), or the AISA Dual airborne sensor (500 bands with spectral resolution of 2.9nm). The performances of traditional serial and parallel algorithms for computing the component trees are strictly correlated to the size and the quantization of the data. The size of remote sensed images is usually in the order of several gigabytes due to the depiction of vast and complex scenes. Consequently, they can not be stored or processed by algorithms designed for a single shared-memory machine. Furthermore, due to the high sensitivity of the new sensors, e.g., radiometric resolution, these images are characterized by an ample domain of integers or floating point values, which directly affect the processing time.

In this paper, a novel shared-and distributed-memory hybrid algorithm for the efficient computation of exact component trees, specifically the min-and max-tree, of integral gray-scale and floating-point images, is presented. For this, the problem, i.e., the image is subdivided into equal-sized chunks that get assigned to all available distributing computing nodes. Then, each of the nodes computes a local, partial component tree of the assigned chunk. A modified version of the shared-memory parallelized, depth-first, flooding maxtree algorithm proposed by Ouzounis et al. [START_REF] Ouzounis | A parallel dual-input max-tree algorithm for shared memory machines[END_REF] is employed. Finally, the obtained partial component trees need to be merged into a correct, monolithic global representation. This is achieved through the iterative resolution and rearrangement of the iso-level edges of the image chunks' boundary trees for each gray-level, marking the major algorithmic challenge. In the proposed approach, the level connections are expressed as tuples-a data structural design that has been used by Flick el al. [START_REF] Flick | A Parallel Connectivity Algorithm for De Bruijn Graphs in Metagenomic Applications[END_REF] for the distributed resolution of genomic graphs. This is a novelty in the mathematical morphology framework and the distributed computation of component trees.

The remainder of this paper is organized as follows. Section 2 provides a brief introduction to component trees. The subsequent Section 3 presents an overview over existing algorithms proposed in the literature. In Section 4 the proposed algorithm for parallel and distributed computation of the component trees is laid out. Complexity considerations and implementation details are explained in Section 5. A study of the algorithms strong and weak scaling as well as comparative study to the current state-of-the-art algorithm is presented in the experimental evaluation in Section 6. Finally, Section 7 concludes the paper, discussing the findings of this work and presents opportunities for future work.

COMPONENT TREES

Component trees were introduced by Jones [START_REF] Jones | Component Trees for Image Filtering and Segmentation[END_REF], [START_REF] Jones | Connected Filtering and Segmentation Using Component Trees[END_REF] as efficient image representations that enable the computation of advanced morphological filters in a simple way. These trees are hierarchical structures that encode the threshold sets and their inclusion relationship. Thereby, each sub-tree is nothing different than all connected image components up until the given gray-threshold. As a result, one of the major advantages of components trees is their possibility to efficiently implement of connected filters.

More formally, let f : Ω → E be a discrete twodimensional gray-scale image, defined on a spatial domain Ω ⊆ Z 2 and taking values on a set of scalar values E ⊆ Z. For any λ ∈ Z, a lower L(f ) and upper U (f ) threshold set is defined by:

L(f ) = {x ∈ Ω, f (x) < λ}, (1)
U (f ) = {x ∈ Ω, f (x) > λ}, (2)
Let P(Ω) be the power set of all the possible subsets of Ω. Given X ∈ Ω, the set of connected components of X is denoted as C(X) ∈ P(Ω). Each connected component is represented by a unique point called the level root [START_REF] Najman | Building the component tree in quasilinear time[END_REF], or canonical element [START_REF] Berger | Effective Component Tree Computation with Application to Pattern Recognition in Astronomical Imaging[END_REF]. Considered two points x, y ∈ Ω, and x r the root of the tree, x is canonical if x = x r or f (parent(x)) < f (x) (where parent is the image that encodes the inclusion relationship of the threshold sets [START_REF] Berger | Effective Component Tree Computation with Application to Pattern Recognition in Astronomical Imaging[END_REF])

If ≤ is a total relation, any two connected components Y, Z ∈ C(L(f)) are either disjointed or nested. The mintree and max-tree structures represent the lower, respectively upper, threshold components in L(f) and U (f) as well as with their inclusion relations. For example, Fig. 1c shows the max-tree structure of the image in Fig. 1a. The arrows denote the parent relation between the nested connected components that are identified in Fig. 1b. This is a simplified case that it is used for clarification purposes. In synthetic images that include more complex shapes (see the example in Section 4) or real scenarios the max-tree structure is less intuitive since its hierarchy is not driven by the inclusion relationship of connected components as it appears in Fig. 1. Furthermore, the notion of tuple as used in this work should be understood as a finite sequence of four elements of the form c x , x, c y , y with x, y ∈ Ω being two pixel coordinates and c x , c y ∈ E ∧ f (x) = c x ∧ f (y) = c y the colors of these pixels.

RELATED WORK

The selection of the most appropriate algorithm for computing the component trees shall be made according to the properties of the input image (i.e., size and pixel value quantization) and the processing resources available such as memory capacity and number of computing cores. Carlinet et al. [START_REF] Carlinet | A Comparative Review of Component Tree Computation Algorithms[END_REF] presented a comparative review of the state-of-theart algorithms and provided detailed guidelines for selecting the most suitable algorithm. The algorithms are grouped into three main classes: immersion-, flooding-and mergebased. Algorithms that belong to the immersion and flooding class, may also be referred to as leaf-to-root merging and root-to-leaf flooding methods, respectively [START_REF] Moschini | A hybrid sharedmemory parallel max-tree algorithm for extreme dynamic-range images[END_REF]. Since this section is not intended to repeat the review, Fig. 2 merely presents a timeline for each algorithm class, and how they have developed in the past years.

As explained in Section 1, the algorithm proposed in this work is of hybrid nature, entailing shared-and distributedmemory parallelization aspects simultaneously. The nodelocal parallelization is thereby based on a flooding strategy, while the distributed computation components follows a merge-based approach. Therefore, a detailed explanation of immersion algorithms is deliberately left out, and the reader is referred to Tarjan [START_REF] Tarjan | Efficiency of a good but not linear set union algorithm[END_REF], Najman et al. [START_REF] Najman | Building the component tree in quasilinear time[END_REF], Berger et al. [START_REF] Berger | Effective Component Tree Computation with Application to Pattern Recognition in Astronomical Imaging[END_REF] and Carlinet et al. [START_REF] Carlinet | A Comparative Review of Component Tree Computation Algorithms[END_REF].

Flooding Algorithms

The first flooding algorithm was proposed by Salembier et al. [START_REF] Salembier | Antiextensive Connected Operators for Image and Sequence Processing[END_REF]. It is an efficient algorithm which retrieves the pixel at the lowest gray-level, i.e., root, through a scanning step and then it performs a propagation by flooding the neighbor at the highest level, i.e., a depth-first traversal of the connected components at higher intensities. Pixels in the propagation front are stored in a hierarchical queue composed by as many First In First Out (FIFO) queues as the number of gray-levels. It allows to directly access any pixel in the FIFO queue at a given level. Salembier's et al. [START_REF] Salembier | Antiextensive Connected Operators for Image and Sequence Processing[END_REF] algorithm was rewritten in a non-recursive implementation by Hesselink et al. [START_REF] Hesselink | Salembier's min-tree algorithm turned into breadth first search[END_REF], later also by Nister et al. [START_REF] Nistér | Linear time maximally stable extremal regions[END_REF] and Wilkinson et al [START_REF] Wilkinson | A fast component-tree algorithm for high dynamic-range images and second generation connectivity[END_REF]. The algorithm presented by Wilkinson aims at solving the limitation of Salembier, the linear scaling with the number of gray-levels. for all p do parent(p) ← -1

⊲ meaning "unseen" ⊲ all points at current level done? q ←TOP(pqueue)

36:

⊲ Attach r to its parent 37:

if f (q) = f (r) then PROCESS-STACK(r, q)

38:

⊲ Particular case of the last element, the tree root:

39:

root ← POP(levroot) Ouzounis et al. [START_REF] Ouzounis | A parallel dual-input max-tree algorithm for shared memory machines[END_REF] 2008 (p) Wilkinson et al. [START_REF] Wilkinson | Concurrent Computation of Attribute Filters on Shared Memory Parallel Machines[END_REF] 2011 (s) Ouzounis et al. [START_REF] Ouzounis | Interactive collection of training samples from the max-tree structure[END_REF] 2017 (p) Moschini et al. [START_REF] Moschini | A hybrid sharedmemory parallel max-tree algorithm for extreme dynamic-range images[END_REF] 2017 (p) Kazemier et al. [START_REF] Kazemier | Connected Morphological Attribute Filters on Distributed Memory Parallel Machines[END_REF] (c) Merge-based Fig. 2: The main three classes of algorithms and their chronology. Each entry includes the publication year, the algorithm class ( s -serial or p -parallel) and the corresponding reference.

a random point start pixel is chosen as the flooding point. This pixel is now considered as a canonical element, i.e., the representative of the connected component, and it is pushed on the stack levroot. This stack stores the representative pixels of the visited components; these pixels are the roots of the sub-trees, so they represent components at different gray-levels. The main purpose of the flooding phase is to compare the gray-level of each pixel p with its neighboring pixels n, and to enqueue those that have not yet been seen. The first processed pixels are p and the canonical element r of its component. These have the highest priority in the queue, i.e., the highest gray-level, and are on top of levroot (p is not removed from the queue). The neighboring pixels n are pushed on the stack only if f (n) > f (p), line 24, which immediately triggers a jump (goto) to the flood label (in line 21). This jump thus emulates a recursive call, which actually corresponds to a depth-first discovery of the tree. At a certain point, all the neighboring pixels of p will be either in the queue or already processed, meaning that the analysis of p has terminated; we cannot progress deeper in the tree. Following this, p is removed from the queue (line 29), parent(p) is set to r, i.e., the canonical element.

In order to ensure that the canonical element will be the last one inserted, p is added to S when r = p (line 31). After p is removed from the queue, the canonical element r is attached to its parent only when the level component has been fully processed (line 30). The last step (startging from line 32) aims at setting the parenthood relationship between components. The first element q of pqueue is retrieved, and the PROCESS-STACK procedure is called (line 37) when q has a different level than p. It pops the stack, sets the parent relationship between the canonical elements, and inserts them in S until the top component has a level no greater than f (q)-lines 3 to 6. When the stack gets empty or the top level is lower than f (q), then q is pushed on the stack as the canonical element of a new component-lines 7 and 8. The top element of the stack is the current root pixel-linelst:last and following. The algorithm ends when all points in queue have been processed, then S only misses the root of the tree, which is the single element that remains on the stack (line 38 and following).

Merge-based Algorithms

The natural way to implement a parallel algorithm is to divide the original image domain and compute the max-tree on each sub-image using any algorithm from Fig. 2. In order to compute this partition, the image should be split in N p connected disjoint regions, which is the union that forms the entire image domain. During this step, the image is split into a reasonable number of chunks, which reflects the underlining processing architecture, e.g., number of threads available. For instance, when the number of image chunks is lower than the number of threads the domain is not decomposed enough and the distribution of the computations is not yet optimal (load imbalance). Some threads will idle while having to wait for other threads to finish. Once all sub-trees are generated, they can be merged into a single global tree as proposed by Matas et al. [START_REF] Matas | Parallel algorithm for concurrent computation of connected component Tree[END_REF], Wilkinson et al. [START_REF] Wilkinson | Concurrent Computation of Attribute Filters on Shared Memory Parallel Machines[END_REF] and Ouzounis et al. [START_REF] Ouzounis | A parallel dual-input max-tree algorithm for shared memory machines[END_REF]. This is a non-trivial phase as it requires that the gray-levels of the connected components are merged and their parent relationships updated. The merging strategy introduced by Matas et al. [START_REF] Matas | Parallel algorithm for concurrent computation of connected component Tree[END_REF] follows the same principle used by Wilkinson. However, the algorithm starts by computing partial 1-D trees (i.e., a tree for each row of the image). Then, the trees that belong to neighboring rows are merged progressively until the global tree is obtained. The merging algorithm proposed by Wilkinson et al. [START_REF] Wilkinson | Concurrent Computation of Attribute Filters on Shared Memory Parallel Machines[END_REF] retrieves the global max-tree and its attribute values from multiple sub-trees that can be derived over any arbitrary image sections. Firstly, each thread computes a data structure (i.e., partial max-tree) that sets the parent pointers and accumulates attribute values. Afterward, all these subdomains are merged through the use of a binary tree. This is achieved through a concurrent merging strategy which connects two partial max-trees step by step. A synchronized mechanism based on two binary semaphores defines when a thread is ready to accept the domain of its neighbor (i.e., the sender has to complete the max-tree computation). Once the last connection is computed by the thread 0, all the threads can resume and proceed with the filtering phase. Ouzounis et al. [START_REF] Ouzounis | A parallel dual-input max-tree algorithm for shared memory machines[END_REF] proposed a max-tree algorithm implementation for attribute filtering based on the concurrent merging strategy in [START_REF] Wilkinson | Concurrent Computation of Attribute Filters on Shared Memory Parallel Machines[END_REF]. The hybrid algorithm proposed in this paper relies on the merging strategy used in [START_REF] Ouzounis | A parallel dual-input max-tree algorithm for shared memory machines[END_REF], [START_REF] Wilkinson | Concurrent Computation of Attribute Filters on Shared Memory Parallel Machines[END_REF] and additional details can be found in Section 4.

Connected Component Labeling

The main problem of generating the max-tree can be split into two parts: find the connected components and establish their hierarchy. The task of grouping connected pixels within an image can be seen as the well-known Connected Component Labeling problem [START_REF] Samet | Efficient component labeling of images of arbitrary dimension represented by linear bintrees[END_REF]. This is an important step for a large number of applications and it relies firstly on finding which parts of an object, e.g., binary images, gray-levels images, data with higher dimensionality, etc., that are physically connected, based on a connectivity rule, and secondly, to label them. Iverson et al. [START_REF] Iverson | Evaluation of connectedcomponent labeling algorithms for distributed-memory systems[END_REF] provide an evaluation of connected-component labeling algorithms in the context of distributed computing, when data that need to be processed in a given application usually require large processing power and distributed-memory machines. They conclude that there is an unavoidable compromise to find between memory and processing time. Most of the available parallel algorithms are problematic especially in terms of memory requirements. For instance, the merging step could end up on a single node resulting in an unbalanced scenario. However, algorithms that try to solve this problem provide poor scaling results in terms of processing time. At the same time Flick el al. [START_REF] Flick | A Parallel Connectivity Algorithm for De Bruijn Graphs in Metagenomic Applications[END_REF] proposed a scalable distributed-memory algorithm to overcome this problems raised by Iverson. They have aimed to solve issues such as excessive memory usage, extra computation and communication of the processors, and load balancing. The main idea of the algorithm is similar to the Shiloach-Vishkin algorithm [START_REF] Shiloach | An o(logn) parallel connectivity algorithm[END_REF] in that it transforms the problem into finding weakly connected components within the Bruijn directed graph [START_REF] De Bruijn | A Combinatorial Problem[END_REF]. It will be shown in Section 4 that the proposed algorithm uses the notion of tuples and inverse doubling in order to connect the overlap zones between the split regions and resolves the connected components and their corresponding parenthood.

DISTRIBUTED COMPONENT TREE ALGORITHM

In the following sections the algorithm for the parallel and distributed computation of the component trees will be introduced. To simplify the explanation, only the max-tree case will be presented. A corresponding min-tree algorithm can be inferred by reversing the order in which the gray-levels are processed. Furthermore, the explanations assume a homogeneous, distributed systems. This refers in particular to the workload distribution of the image into contiguous, equalsized partitions. For heterogeneous, distributed systems a different, more appropriate strategy must be chosen, such as proposed by Qin et al. [START_REF] Qin | A Dynamic and Reliability-driven Scheduling Algorithm for Parallel Real-time Jobs Executing on Heterogeneous Clusters[END_REF].

Definitions and Notation

A two-dimensional gray-scale image f can be seen as an undirected graph G = (V, E). V represents a set of verticesthe pixels of the image-and n = |V | the total number of pixels. Then E, a number of edges or non-ordered pairs of vertices (v i , v j ), with i, j ∈ [0, n[, which model the neighborhood relationship of the pixels. Classically, images are either four-or eight-connected [START_REF] Soille | Morphological Image Analysis: Principles and Applications[END_REF], meaning the top, left, right and bottom neighbors, respectively including the diagonals, are considered connected neighbors. The entire graph G is said to be connected if, for any p, q ∈ V , there exists a path from p to q, which is a sequence of s > 1 vertices-i.e., p = p 1 , ..., p s = q-such that every p i ∈ V , and any two successive pixels of the sequence are adjacent e pi,p i+1 ∈ E. Given this definition a connected component

CC is a subgraph of G such that V CC ⊆ V G , E CC ⊆ E G , ∀ p ∈ V CC : f (p) = c, that is maximal: ∀ e = (p, p ′ ) ∈ E CC such that p ∈ V CC and p ′ ∈ V CC , f (p ′ ) = c.
A connected component can be either weak or strong connected, depending on the path length s. Weak connected graphs can have an arbitrary path length, while for strong connected graphs s = 2 holds. Furthermore, if not stated otherwise, the following symbols are defined for the remainder of the document: h and w is the height and the width of the image f , respectively. The entire image has a gray-level depth d, i.e., the number of different gray-values c. Concerning parallelization, the number of available distributed compute nodes is p, while the local number of shared-memory threads is labeled with t. For the explanation of the distributed resolution, it is also necessary to introduce what is coined a tuples. These are essentially quartuples, mathematical tuples with four components, of the form c i , p k , c j , p l with c i and c j being two gray-values and p k and p l two vertices. They are used to explicitly express edges e ∈ E of the image f of a canonical point with a certain gray-value to a different canonical point of a given other gray-value.

Concept

Algorithm 2 Pseudo-code of the proposed distributed maxtree algorithm. t ← number of threads 6: tuples ← RESOLVE(area tuples, root tuples)

f ′ ← LOAD-PARTIAL-IMAGE(f, r, p

return parents

The general nature of the distributed max-tree algorithm can be described as divide-and-conquer. This means, that the entire problem, i.e., the image, is divided into sub-images for which the respective max-trees are computed and that are then successively merged along the division boundaries. The major algorithmic challenge lies in the latter stage. It requires to solve two demanding graph theory sub-problemsconnected component labeling and graph canonicalizationin distributed memory environments. This work proposes an iterative, parallel merging algorithm based on explicit expression of the max-tree edges as directed tuples. For this, halo-zones are employed-one-pixel wide, redundant overlaps of the partial images of neighboring image chunks (see also Figure 3). Conflicts within the tuples signify the need to rearrange the edges of the boundary tree in order to obtain a globally correct view. Such a conflict could for example be, that the exact same gray-level component points to different parent components in the partial trees.

Analogous to the two graph problems to be solved, the proposed algorithm requires two kinds of tuples. First, the locally determinable edges of the boundary trees in the halo zone (root tuples), and, second, information about components that have been split due to the image division (area tuples) and their the canonical points. In the proposed resolution approach, the tuples are iteratively scanned for and remapped to the most optimal candidates. This entails merging split iso-level and determining the best parent or root each component and as a result remapping the tuples. Less optimal tuple candidates are replaced by transient edges to ensure a correct merge of the remaining trees. Subsequently, the computed changes need to be applied to the partial images. For this, the resolved tuples are send back to the respective sub-images of origin, and utilized to obtain the globally create max-tree. Algorithm 2 sketches the proposed strategy.

Local Max-Tree Algorithm

For the local computation one can in principle employ any correct max-tree algorithm. The proposed solution specifically utilizes a modified version of the recursive Salembier's depth-first flooding algorithm [START_REF] Salembier | Antiextensive Connected Operators for Image and Sequence Processing[END_REF]. As observed by Carlinet et al. [START_REF] Carlinet | A Comparative Review of Component Tree Computation Algorithms[END_REF], Salembier's algorithm was rewritten in a nonrecursive implementation by Hesselink et al. [START_REF] Hesselink | Salembier's min-tree algorithm turned into breadth first search[END_REF] and later by Nister et al. [START_REF] Nistér | Linear time maximally stable extremal regions[END_REF] and Wilkinson et al. [START_REF] Wilkinson | A fast component-tree algorithm for high dynamic-range images and second generation connectivity[END_REF]. In [START_REF] Wilkinson | A fast component-tree algorithm for high dynamic-range images and second generation connectivity[END_REF], it was shown that replacing the hierarchical queue (used in [START_REF] Salembier | Antiextensive Connected Operators for Image and Sequence Processing[END_REF]) by a priority queue to perform the recursive flooding strongly reduced the computation of component trees, especially for high-dynamic range images (i.e., with floating point). Due to the similarities of approach and the proved efficiency in terms of processing time in [START_REF] Wilkinson | A fast component-tree algorithm for high dynamic-range images and second generation connectivity[END_REF] and [START_REF] Nistér | Linear time maximally stable extremal regions[END_REF], Carlinet et al. [START_REF] Carlinet | A Comparative Review of Component Tree Computation Algorithms[END_REF] merged these solutions and suggested a novel non-recursive implementation, which is utilized in this work. Additionally, it has been enhanced to always use the minimal pixel index as canonical point for an iso-level and to yield better computational performance. Algorithm 3 presents the corresponding pseudo-code.

The first change can best be seen in line [START_REF] Aldeborgh | Unsupervised object detection on remote sensing imagery using hierarchical image representations and deep learning[END_REF]. In contrast to the original non-recursive variant, the canonical area point is not chosen at the beginning of an iso-level processingas it may not yet be the canonical minimum-but rather constantly maintained throughout the process. This is done by keeping the current area minimum at a specific place, e.g., the front of the pixel vector, and compared to on insertion of new elements. Only after all pixels of the entire iso-level is found, the canonical point is assigned in the parent image, see also line 31, and thus minimality of the index guaranteed.

Moreover, when one considers the computational performance of the algorithm, the proposed modifications also allows for faster computation. Before, each gray-level had its one hierarchical queue in Salembier's original algorithm formulation or a singular in Carlinet's non-recursive reformulation. This approach scales logarithmically with both, the number of gray-levels as well as the number of pixels per channel. In the proposed variant the gray-levels are keys to a map, called stacks, that has vectors for the corresponding pixels as keys. Then, insertions only scale logarithmically with the number of gray-levels, for locating the respective vector in the map, but the actual push operation happens in constant Algorithm 3 Pseudo-code of the modified version of Salembier's depth-first, flooding-based max-tree algorithm. The strategy for shared-memory parallelization explained by Ouzounis et al. [START_REF] Ouzounis | A parallel dual-input max-tree algorithm for shared memory machines[END_REF] has been chosen. Their max-tree algorithm was based on the parallel implementation proposed by Wilkinson et al. [START_REF] Wilkinson | Concurrent Computation of Attribute Filters on Shared Memory Parallel Machines[END_REF]. The authors largely proved the correctness and efficiency of the approach and they predicted an achievable (near) linear speed-up beyond 4 CPUs. In greater detail, the local image partition of the node distribution step, is again virtually partitioned in t equally sized horizontal chunks, without overlap, and assigned to one of the t threads. For each of the partitions the partial max-tree is computed using the introduced algorithm. The virtual images boundaries are realized by excluding respective pixels in the neighborhood searches. Finally, the partial max-trees are merged using the connect function, equally presented by Ouzounis et al [START_REF] Ouzounis | A parallel dual-input max-tree algorithm for shared memory machines[END_REF]. Minor changes have been made to ensure that the canonical points of each iso-level is guaranteed to be minimal. This can be achieved by performing a look-ahead on the upcoming elements of the merge stacks and potentially swap them if required. After the local computation, one would obtain the partial max-trees depicted in Fig. 3e and Fig. 3f.

Tuple Generation

The tuple generation of the partial images is subdivided into two major steps. First, there is the generation of root tuples, representing the branches of the halo areas' partial maxtrees, and, second, the area tuples-i.e., tuples that connect or stitch the gray-level components divided by the image partition.

The former set of tuples is generated by traversing the boundary max-trees recursively upwards towards the root. For each edge within the tree a corresponding tuple is created with the canonical point of a sub-ordered, child iso-level pointing to the canonical point of its parent. In line with the used notation it is the result of a HALO-TREE-EDGES function invocation. The performance of this step can be further optimized by skipping already visited branches also avoiding redundant tuples. Generating the root tuples can be performed fully local and does not require any data exchange with other processing nodes.

The former step, the connection of the split components, is achieved by identifying differences in the labeling of the canonical points of said components in the halo areas. It is the result of a HALO-COMPONENTS function invocation. This can achieved efficiently by performing two prefix-sums pixels in the halo area of the partial images. In the first run, the canonical points are broadcast downwards across the partial images, while the second prefix-sum, in reverse, back-propagates found differences. Each conflict implies that a merge of the two components is necessary. Therefore, an area tuple needs to be created. It maps the higher canonical point x h to the lower x l , resulting in a tuple of the form c, x h , c, x l . For a faster collapse of transitive chains in the resolution stage, the inverse tuple c, x l , c, x h is additionally generated and stored. This entire process of performing the prefix-sums can be efficiently realized using a logarithmic merge tree across all available nodes.

Distributed Tuple Resolution

The resolution of the tuples happens essentially in the same mode as the generation of tuples. First, the area tuples are resolved, which requires the iterative resolution of the weakly connected area components into strong ones, and then, second, the resolution of the edges for the normalized components. The respective high-level pseudo-code is displayed in Algorithm 4. It makes use of the distributed communication primitive ALLREDUCE several times. In line with this text it should be understood as a function that reduces, i.e. combines, vector of values with the same length and arbitrary, but consistent data types, element-wise using an operator given as second parameter. The single, final result vector is broadcast to all participating distributed compute nodes. ALLREDUCE can be efficiently implemented using a logarithmic communications tree.

Algorithm 4 Pseudo-code of a single iterations of the distributed area tuple resolution. GLOBAL-SORT(area)

13:

rules ← RESOLVE-COMPONENTS(area(grayv))

14:

unresolved ← REMAP(area(grayv), rules)

15:

unresolved ← ALLREDUCE(unresolved, Or)

return tuples

In general, the area tuple resolution is inspired by the distributed connected component labeling algorithm presented by Flick et al [START_REF] Flick | A Parallel Connectivity Algorithm for De Bruijn Graphs in Metagenomic Applications[END_REF]. The goal is to turn a weakly connected graph, here the components, into a strongly connected ones, meaning directly pointing to the correct canonical point of an area without intermediate, transitive connections. This is achieved by following the graph edges until the most optimal, i.e., smallest pixel point is found. However, in distributed memory environments it might not be possible to follow all paths directly as they might be residing on a different machine. To overcome this, the longest local partial graph paths are computed and iteratively shortened until the strong graphs are found. For this, all tuples are globally sorted, i.e., across all nodes, and locally linearly scanned for the most optimal candidate, remapped and saved. This process is repeated until convergence is achieved, which is equal to having no tuples remapped in the current iteration. Technically, there are two challenges involved in this.

First, there is the problem of globally sorting (GLOBAL-SORT in Alg. 4) the tuples. This means that all tuples need to be partially ordered, so that the smallest element is on the node with the smallest rank and the maximal tuple on the node with the highest rank. For this, the distributed max-tree algorithm uses an enhanced version of the parallel sorting by regular sampling algorithm [START_REF] Shi | Parallel Sorting by Regular Sampling[END_REF]. In the variant that is proposed here, the number of tuples are additionally balanced after sorting them, in order to keep the workload equal on each node. Algorithm 5 Pseudo-code of the component resolution algorithm mapping weakly connected tuples into stronger ones.

return rules

Second, there is the problem of finding the strongest connected components local to each node. The approach proposed here, is sketched in Algorithm 5. Each of the area tuples is linearly scanned and the most optimal, i.e. smallest canonical point memorized in an associate map called rules. These can then be applied to the tuples by scanning them once again and modifying the During the transitive solution the entire rule chain needs to always be canonized and the currently known minimum pointed to (lines 8-12). If the tuples are balanced during the global sorting step, the resolution of the partial component graphs, i.e., all the tuples with the same origin, may additionally be fragmented across the memory of multiple nodes. Therefore, the distribute max-tree algorithm must connect these partial graphs in each iteration. This can be achieved by exchanging the start and end of the sorted tuple chain including the found canonical point with the direct neighbors. Given that the neighboring node proposes a better canonical point, it is adopted instead of the one found locally. Transitivity is achieved by logarithmically merging these chains across all available machines. For this two prefix-sums (sometime also cumulative sum or scan, lines [START_REF] Benediktsson | Hierarchical Analysis of Remote Sensing Data: Morphological Attribute Profiles and Binary Partition Trees[END_REF][START_REF] Berger | Effective Component Tree Computation with Application to Pattern Recognition in Astronomical Imaging[END_REF], first from left to right across the ranks, and then in reverse to propagate potentially better canonical points back.

Furthermore, each tuple exists twice. Once in the "correct direction" pointing from the larger pixel index to the lower index and its inverse, pointing from low to high. Whenever a tuple is remapped, the tuple is flipped, i.e., the direction is changed, for the next iteration in order to back propagate this change to its inverse. The reason behind this is, that the inverse tuple might have found an even more optimal canonical point coming from the other side of the chain, due to say half circular structures on the image. Both tuples are then updated and the transitive chain collapsed much quicker. As an effect of this, the number of iterations heavily reduces, as introduced by Flick et al [START_REF] Flick | A Parallel Connectivity Algorithm for De Bruijn Graphs in Metagenomic Applications[END_REF].

After the canonicalization of the weak connected components into strong connected components, the roots for each of the canonical points must be found. The corresponding approach is sketched in Algorithm 6. For this, the area tuples and root tuples are merged first in order to normalize of the components' canonical points (lines 2-4). Only then can the most optimal root for each component be determined. This is achieved by once again linearly scanning the tuples and memorizing the best candidates in an associative map called best roots. Similarly to the area tuples fragmentation, root candidates of a single component may be scattered across multiple nodes. Therefore, they must be connected two prefix-sum operations (lines 7-10).

In a third and final linear scan over the combined tuples the found best roots are evaluated. There are four possible options. First, the tuple's root gray-value is smaller, i.e. it is further up in the tree, than the best root. In this case, a tuple needs to be created that connects the best root transitively with the one from the tuple (see case 2). Second, the grayvalue of the tuple and the best root match, but the root has a smaller canonical point (see case 3) This means that the root and its connected component as well as the one of the neighbor are weakly connected via the area of the current tuple. An area remapping tuples needs to be created, including its inverse, and added to the area tuples. Third, the current tuple already points correctly to the best root (see case 4), then the tuple is inverted, essentially pointing the "wrong" way around from the root to the lower area and pushed into the root's gray-valued tuple bucket. The reason behind this is, that the canonical point of the root might still be changed during the area resolution phase of its gray-value. Therefore, it may not be marked as finished yet, but kept until the respective gray-level of the roots has been resolved. Finally, the fourth condition (see case 1) is meant for inverted root tuples. After their normalization, they are flipped yet again into the "correct" order, i.e., pointing from high grayvalues to low gray-values and send back to its respective tuple bucket.

Obtaining the Global Parent Image

After the tuples have been resolved the global parent image can be obtained by redistributing the tuples back to their original sub-image. For this, the each tuple is send back to the node, that contains the pixel with the index of the second tuple component, independent whether it is a area or edge tuple. Each of the area tuples is then stored in an associative container, mapping from the original points to the canonical point. Subsequently, while iterating over the image, each of the pixels is normalized to its correct canonical points using said data structure. The root tuples are handled slightly differently. Each tuples is visited once and the pixel index at the from part of the tuples is simply set to the destination of the tuple. After this step, each node posses the correct partial parent image representing the global max-tree. 3) Iteratively resolve the generated the tuples beginning from the highest to the lowest gray-level-component first, then edges. Tuples with ↓ are dropped (incorrect or inverted), mark conflicts during remapping and applied changes are bold.

RESOLVE C3 RESOLVE-COMPONENTS NORMALIZE RESOLVE-ROOTS Generated tuples [ 3, 3, 3, 5 ↓, 3, 3, 3, 9 ↓] [ 3, 3, 0, 0 , 3, 5, 1, 1 , 3, 5, 3, 3 ] [ 3, 3, 0, 0 , 3, 3, 1, 1 ] [ 2, 19, 0, 0 , 2, 19, 1, 1 ] [ 3, 5, 3, 3 , 3, 9, 3, 3 ] [ 3, 9, 2, 19 , 3, 9, 3, 3 ] [ 3, 3, 2, 19 ] [] C2 RESOLVE-COMPONENTS NORMALIZE RESOLVE-ROOTS Generated tuples [] [ 2, 19, 0, 0 , 2, 19, 1, 1 ] [ 2, 19, 0, 0 , 2, 19, 1, 1 ] [ 1, 1, 0, 0 , 1, 1, 1, 18 ] [] [ 2, 19, 1, 18 ] [ 2, 19, 1, 18 ] [ 1, 18, 1, 1 ] C1 RESOLVE-COMPONENTS NORMALIZE RESOLVE-ROOTS Generated tuples [ 1, 1, 1, 18 ↓] [ 1, 1, 1, 18 ] [ 1, 1, 0, 0 ] [ 0, 0, 0, 8 ] [ 1, 18, 1, 1 ] [ 1, 18, 0, 8 , 1, 18, 1, 1 ] [ 1, 1, 0, 8 ] [ 0, 8, 0, 0 ] C0 RESOLVE-COMPONENTS NORMALIZE RESOLVE-ROOTS Generated tuples [ 0, 0, 0, 8 ↓, 0, 0, 0, 8 ↓] [ 0, 8, 0, 0 ] [] [] [ 0, 8, 0, 0 , 0, 8, 0, 0 ↓] [] [] []
Fig. 3: Toy example demonstrating the tuple resolution of the distributed max-tree algorithm.

(k × d × w d × log(w × d) × log(p)) O(w × d) O(k × d × log(p)) Redistribution O(w × d) O(w × d) O(log(p)) Application O(w × d + n p × log(w × d))) O(w × d) -
Algorithm 6 Pseudo-code of the distributed resolution of the root tuples.

1: function RESOLVE-ROOTS(grayv, area, roots)

2:
combined ← CONCAT(area(grayv), roots(grayv))

3:

GLOBAL-SORT(combined) LEFT-PREFIX-SUM(ends, Min) return tuples

IMPLEMENTATION

The proposed parallel algorithm has been implemented in C++ and is available on the open-source code repository Github [START_REF] Goetz | Distributed Max Tree Implementation[END_REF]. The coarse-grained parallelization across multiple nodes has been realized using the Message Passing Interface (MPI) [START_REF] Gropp | Using MPI: Portable parallel programming with the message-passing interface[END_REF]. For the shared-memory implementation C++11 native threads have been used. The algorithm accepts data loaded from files in the Hierarchical Data Format 5 (HDF5) format [START_REF]Hierachical Data Format 5[END_REF], which it will also store the resulting parent image to.

Complexity

In this section the time and space complexity for the algorithm steps of the distributed max-tree computation are laid out. A summary can be found in Table 1. The used symbols are explained in Section 4.1. All formulas are given for the worst-case scenario.

The time and space complexity for loading the sub images can be straight-forward inferred and amount to the number of total pixels divided by the amount of available processing nodes, as each of the receives an equally size chunk of the entire problem. This requires to potentially exchange one message in which the image dimensions are broadcast. For the local max-tree computation each of the t thread needs to allocated the part of parent image, that is equal in size to the processed raw image, plus and additional area remapping that is solely dependent on the image width, explaing the space complexity. The computational complexity consists of the linear image-scan for each thread and sub-image, which in turn need to do look-ups into the associative container for the stacks, resulting in the first summand. However, each thread needs to be merged with its direct neighbors, which can be done in logarithmic fashion as explained in Section 4.3, along the virtual split boundaries, i.e., the width of the image.

The next two steps involve the generation of the tuples. In the worst case, for each of the boundary pixels-again the width of the image-a tuple needs to created, resulting in according space and time complexity. Although, in practice the number will most of the times be much lower, because the boundary zones mostly consists of connected flat zones with shared parents nodes, resulting in early outs. In fact, the average complexity should therefore be closer to O(w * log( d 2 )), but is for obvious reasons dependent on the analyzed data. Contrary to root tuples, which can be generated entirely local, area tuples needs to be stitched together across all processing nodes, resulting in the additional factor for the time complexity. Using a two prefix-sum operations this information can be exchanged in a logarithmic number of communication steps, explaining the messaging complexity.

The main resolution consists of k iterations of sorting, linearly scanning and remapping the tuples for each of the available gray-levels of the gray depth d. It is assumed here that the number of tuples per gray-level is more or less evenly distributed, resulting in the purposefully chosen term w d (it would actually cancel out with d). The sorting adds both logarithmic components, over the number of tuples and nodes, as it requires reordering them across all machines. One of the major uncertainty factors is the iteration constant k, which is dependent on the data. In the worst case, k is equal to the number of processing nodes p as upper bound, given a single tuple needs to visit each single machine due to transitivity. However, in practice, the iteration count will remain low, typically only one or two, even for a high number of nodes, due to the area stitching step and tuple inversions, effectively minimizing the canonicalizations. The corresponding message complexity can be explained in a similar fashion. For each of the d gray-levels, k iterations are performed requiring the communication of the tuples using a logarithmic communication primitive across all cores. Thereby, the tuples can be exchanged in whole, not requiring to break them down into individual messages.

Finally, the resolved tuples need to be send back to the partial image of origin and locally applied. The former step, requires to exchange messages between each of the nodes. Using a logarithmic communication tree, this can be done ad hoc, explaining the message complexity. The later step requires iterating over the whole image and normalizing each of the pixels using an associative data structure with logarithmic look up time. This is the reasons for the second summand in the time complexity equation. The first can be explained by the changing roots by directly assigning them while iterating through the tuples.

Generally, the algorithm has complexity classes that are either linear or linear-logarithmic, supporting good scalability overall. The only bottleneck seems to be the iterative constant k that could potentially degrade into the number of used compute nodes p as upper bound. However, this is rarely the case in practical use and presents an opportunity for future research.

Implementation Details

For the algorithm implementation the Message Passing Interface (MPI) [START_REF] Gropp | A high-performance, portable implementation of the mpi message passing interface standard[END_REF] programming framework has been used. It provides low-level network communication primitives to exchange one-to-one and many-to-many messages between the participating distributed nodes. Efficient algorithms usually rely on the later category of operations, the so-called collectives, due to possibility of achieving a logarithmic scaling, in time and number of messages, across the number of cores. The proposed implementation makes use of two concrete function. On the one hand this is MPI_Allreduce and on the other hand this is MPI Scan. The latter is a concrete realization of the prefix-sum operation, accepting a vector of input values and element-wise calculating the left partial sum with respect to the passed binary operator. In the proposed algorithm implementation is used at various points, e.g. connecting the components in the halo areas or to exchange the canonical points of the partial strong connected components (see Section 4.5).

Specifically for the generation of the area tuples that connects the halo zones, two prefix-sums, alternating from left to right, are necessary. The first operations identifies conflicts in the canonical point labeling and the second propagates them back. In order to retain the complete information about all remapping rules for all nodes, the prefixsum function (MPI_Scan) would require an exchange buffer with a worst case memory complexity of O(p * w * 2)-i.e. two halo zones 2 * w for each of the p nodes. This is highly undesirable, as it scales both, with the number of processing nodes as well as the width of the image. One can realize this operation more efficient, if only the outer boundaries of the already merged images are communicated and the intermediary remapping rules are memorized in a data structure, e.g., a map, across the two MPI_Scan calls. Then, the memory complexity of the sent buffer simply becomes O(w * 2).

The MPI framework does allow the registration of custom functions for collective calls, such as MPI_Scan. These must be associative and optionally commutative, which is satisfied by the above operation. In practice, however, such a reduction function additionally needs to have static linkage, or in other words, it must be a singleton. Furthermore, due to the definition of the MPI API standard, it is also not possible to pass any context or state, say an object pointer to the aforementioned map, to the prefix-sum. For this reason, a straight forward realization of the stateful MPI_Scan is not possible. Yet, there is the alternative of working around this limitation by accessing static data structures, e.g. a map modeling the local execution context, within the scope of the custom reduction operation. This could be a globally defined variable or static class member.

Algorithm 7

Pseudo-code of a thread-safe, stateful MPI reduction operation and subsequent usage by a prefix-sum. In this case, though, the whole distributed max-tree implementation effectively also becomes a singleton and may not be used in multi-threaded environments, which is suboptimal for a number of analysis uses cases. Therefore, it has been chosen an approach as sketched in Algorithm 7. A set of potential remapping data-structures from different threads is stored in a global associative container, here rules, guarded by a mutex. Before calling an MPI_Scan the remapping data-structure must be stored and during the custom reduction operation retrieved. In order to be able to correctly retrieve the remapping data-structure a unique, shared key must be chosen, for example the current thread identifier.

One enhancement to and possibility for future research on the MPI standard versions could be, to directly allow passing context pointers to every API call that utilizes reduction operations. This pointer is simply forwarded to the custom reduction operation on each invocation and then used to realized stateful behaviour. In case a context is not needed, it may be set to null. As a result, this would remove the locking overhead and the code becomes cleaner and more understandable.

EXPERIMENTAL EVALUATION

Environment

The experiments have been performed on the JURECA system [START_REF]JURECA: General-purpose supercomputer at Jülich Supercomputing Centre[END_REF] at the Juelich Supercomputing Centre. The system consists of 1884 compute nodes with each having two Intel R Xeon R E5-2680 v3 Haswell CPUs with 12 cores at 2,5 GHz and Hyperthreading. 1604 compute nodes have 128 GiB, 128 nodes 256 GiB, 74 node 512 GiB and two nodes 1024 GiB DDR5 RAM. For our experimental evaluation the following software libraries have been used-HDF5 1.8.18 parallel and ParaStation MPI 5.1.9. Source code has been compiled with g++ 5.4.0 optimization level 03. The available benchmark for the experiments relies on a maximum of 32 nodes and 24 threads.

Datasets

As for the used data, the tests have been performed on two real-world images depicted in Figure 4. The first dataset is a Pléiades Ortho Product 1 that was acquired over the Naples metropolitan area (Italy) in 2013. It includes four Pansharpened images with spatial and radiometric resolution of 0.5m and 8 bpp, respectively. This dataset was selected due to its free availability, its sufficient size and spatial resolution, which are relevant to the needs of the remote sensing community for scalable and accurate methods that allow to classify rapidly and accurately objects of interest over vast areas, as was discussed in Section 1 (e.g., produce very high resolution land cover mapping at the European scale). The second dataset is an image that was taken at the ESO Paranal Observatory in Chile by the Visible and Infrared Survey Telescope for Astronomy (VISTA). It portraits more than 84 million stars in the central regions of the Milky Way [START_REF] Saito | Milky Way demographics with the VVV survey[END_REF]. The Figure 4a and Figure 4b show the true-color image of both datasets. For the Naples dataset, experiments are performed only using the first channel [START_REF]Pleiades Ortho Product Dataset: Pan-Sharpened Green Band Channel[END_REF]. For the ESO, the RGB image is simplified to a singular luminance channel, similarly to how it was done by Moschini et al. [START_REF] Moschini | A hybrid sharedmemory parallel max-tree algorithm for extreme dynamic-range images[END_REF] in order to obtain similar conditions for benchmarking the proposed algorithm. The luminance image is obtained through weighing and summing the channels, so that L = 0.2126R + 0.7152G + 0.0722B. However, in order to show that the algorithm scales regardless of the domain size of the gray-levels, three different quantization levels are derived from the luminance channel: 8-uint bpp, 16-uint bpp and 32-float bpp [START_REF]ESO Dataset: 3 Luminance Channels with Different Quantization Levels[END_REF]. Contrary to [START_REF] Moschini | A hybrid sharedmemory parallel max-tree algorithm for extreme dynamic-range images[END_REF], the original size of the image is preserved (≈9Gpx), since the JURECA system provides node with large memory.

Experimental Setup

As discussed in Section 3 there are a number of other serial and parallel versions of the algorithm. Most of them report 1. http://www.intelligence-airbusds.com/en/23-sample-imagery different value permutations for the computation time, memory consumption, speed-up and scalability of their implementations. Carlinet et al. [START_REF] Carlinet | A Comparative Review of Component Tree Computation Algorithms[END_REF] provide their used benchmarks, datasets and the source codes in C++ for many different serial and parallel algorithms 2 . In order to compare results achievable by using serial and parallel computing, the Berger et al. [START_REF] Berger | Effective Component Tree Computation with Application to Pattern Recognition in Astronomical Imaging[END_REF] algorithm has been selected. Moschini et al. [START_REF] Moschini | A hybrid sharedmemory parallel max-tree algorithm for extreme dynamic-range images[END_REF] proved that Berger is the fastest sequential algorithm for images with high quantization values and even floating. However the algorithms depend on the MILENA image processing library [START_REF] Levillain | Why and how to design a generic and efficient image processing framework: The case of the Milena library[END_REF] (i.e., provide fundamental image types and I/O functionality) which was not designed to handle very large images and floating values. For these reasons it was necessary to re-write a new C++ implementation of the algorithm which is library independent. For the parallel processing case, the hybrid shared-memory parallel max-tree algorithm developed by Moschini et al. [START_REF] Moschini | A hybrid sharedmemory parallel max-tree algorithm for extreme dynamic-range images[END_REF] was considered. The algorithm has been implemented in C using POSIX threads and the source code is available publicly 3 . Contrary to the MILENA library, this algorithm has been proposed with the purpose of dealing with large-scale and high-dynamic range images, and was therefore ready to be used out-of-the-box. It may be argued that the comparison is not entirely fair due to the different nature of the algorithms-i.e., shared-memory and distributed-memory-but can very well be investigated for the same number of utilized cores. The expectation naturally is that distributed memory implementation, as the one proposed here, are naturally going to have more overhead compared to shared-memory versions. To the best of our knowledge the only distributed max-tree algorithm has been proposed recently by Kazemier et al. [START_REF] Kazemier | Connected Morphological Attribute Filters on Distributed Memory Parallel Machines[END_REF], but the source code was not obtainable at the time of writing as it is not yet released. The performance assessment of the algorithms proposed by Berger and Moschini against the algorithm proposed in this paper is conducted with two kinds of benchmarks. The first type is focused on the computation time and speedup, while the second measures memory consumption. Each benchmark configuration, meaning a particular node and core count, is executed five times and the following statistics are reported: mean µ, standard deviation σ, minimum, maximum, and coefficient of variation (CV ), defined as ν = σ µ [START_REF] Abdi | Coefficient of variation[END_REF]. The use of the multithreading/MPI hybrid features of the algorithm allows to span the MPI process on each node available and to parallelize it locally using multithreading. For this reason, both types of benchmarks are performed on each number of cores, as shown in Table 2. The strategy is to evaluate first the performance on one core of one node. Afterwards the number of threads are doubled alternating with doubling the number of nodes, until the maximum of 128 cores across nodes and threads is reached.

Speed-up and Memory Consumption

The Figure 5 depicts the experimental results related to the processing time. For each dataset, the plot of the mean execution time and the plot of the speed-up for increasing number of cores is reported. In order to make a fair comparison with the state-of-the-art results (i.e., the shared-memory algorithm [START_REF] Moschini | A hybrid sharedmemory parallel max-tree algorithm for extreme dynamic-range images[END_REF]), the proposed algorithm is first run on a single node. The algorithm's execution time measures the beginning and end of the main() function of the process with the MPI rank 0 and the thread number 0. The speedup coefficient is computed as t p = t 1 /t c . the fraction of the execution time with a single core and the execution time with multiple processing cores. Generally it can be said that the proposed algorithm is able to gain a substantial speed-up for both data sets and the different gray-levels quantizations.

For the 8bpp case, the algorithm shows a constant, near linear speed-up curve. In both datasets, the speed-up shows an increasing behavior for up to 256 cores with no reason to doubt its consistency for a higher number of cores, with an execution time of 9.38 and 62.35 seconds for Naples and ESO, respectively. However, for ESO 16bpp and 32bpp the speed-up flattens sooner, stabilizing at 64 and 16 cores, respectively. With a high gray-level depth, the number of tuples is increasing sharply, resulting in larger merge time. The effect observable here is the Amdahl speed-up boundary for a constant workload. When these results are compared with the Moschini algorithm, the proposed algorithm always provides faster execution times. Unfortunately, for the dataset ESO 32bpp it was not possible to derive any conclusions since the Moschini algorithm did not terminate. A more detailed analysis of the execution time for the different phases of the algorithm (see the Algorihtm 2) is depicted in Figure 6. The results are related only to a single dataset case (ESO 16bpp) because of space considerations. Each set of rows depicts a a specific number of nodes (i.e., 1 node, 2 nodes, 4 nodes and 8 nodes). The time distribution for increasing number of threads is shown in each row. For the single node case, the computation of the local max tree is the most time-consuming phase. This is a shared-memory scenario where the three phases concerning the management of the tuples do not take place. When the number of threads increases, the second most costly phase is the local apply. The local merge needs to be considered only beyond four threads. The same conclusions can be derived for the remaining nodes configuration. However, since it is a distributed memory environment, the phases connected with tuple handling are also present.

In the two-node case the tuple generation phase and the global apply are mostly present, the weight of the tuple resolution becomes more pronounced the higher the node count. This behavior can best be explained by the algorithm complexity, explained in Table 1, showing a logarithmic scaling with the number of processing cores. As has been shown here, the parallel implementation allows to achieve a significant processing time gain when compared with serial processing. In Table 3 the processing times for the different datasets are presented. When considering ESO 32bpp, which is the more challenging dataset used in this work, the proposed algorithm computes the max-tree in ≈27 minutes (with 24 cores) while Berger converges only after ≈3 hours.

Last but not least, comments should be made regarding the memory consumption of the proposed algorithm. Considering the usual trade-off between memory consumption and computational time, the experiments show the proposed algorithm is more memory efficient and takes shorter computational time than Berger and Moschini. Table 4 and Table 5 scrutinize the memory consumption (in GB) for the different algorithms with Naples and ESO 16bpp images, respectively, when computing on a single node. For each given number of threads, it can be noticed that the average and the maximum memory usage of all the tasks in the job are always lower for the proposed algorithm. This means that the algorithm is able to scale in terms of memory consumption and communication cost with respect to large datasets and the number of parallel cores. This is an important factor, considering most of the time the main constraint lies in the memory size.

CONCLUSION

In this work a new parallel and distributed algorithm for the computation of the max-tree of an image has been presented. The parallelization strategy consists of splitting the entire problem, i.e., the image, into equal-sized sub images, for which the partial max-trees are computed that are subsequently merged at the split boundaries. Using this algorithm, substantial speed-ups and scalability could be achieved in computing the max-tree on large real-world images, outperforming the state-of-the-art shared memory implementation. In particular, faster execution time and significantly less memory consumption can be achieved. The proposed algorithm allows to process gray-scale image of arbitrary gray-level depth including floating point values. This makes it suitable for the usage in large-scale image classification task, such as land cover type prediction, which is one of the major practical application domains.

In future work, the equivalent min-tree algorithm including distributed attribute filter are going to be implemented. This will set a solid foundation for the next research goal, the massive parallelization of the tree of shapes [START_REF] Crozet | A first parallel algorithm to compute the morphological tree of shapes of nD images[END_REF]-a contrast independent component tree representation of images.

1 : 7 :

 17 procedure PROCESS-STACK(r, q) while levroot not empty and λ<f (TOP(levroot)) do 5: INSERT FRONT(S, r) 6: r ← parent(r) ←POP(levroot) if levroot empty or f (TOP(levroot)) = λ then

1: @parallel 2 :

 2 function DISTRIBUTED-MAX-TREE(f ) 3: p ← number of nodes 4: r ← processor id in range [0, p[ 5:

  )

7 : 8 :

 78 parents ′ ← LOCAL-MAX-TREE(f ′ , t) 9: 10: root tuples ← HALO-TREE-EDGES(f ′ , parents ′ ) 11: area tuples ← HALO-COMPONENTS(f ′ , parents ′ ) 12:

1 : 7 :

 17 function RESOLVE-COMPONENT(area tuples) if to > f rom then SWAP(to, f rom) 8: canonical ← CANONIZE(rules, f rom) 9: min ← MIN(canonical, to) 10: max ← MAX(canonical, to) 11: rules[f rom] = min 12: rules[max] = min 13: 14: ends ← [FRONT(area tuples), BACK(area tuples)] 15: LEFT-PREFIX-SUM(ends, rules, Min) 16:RIGHT-PREFIX-SUM(REVERSE(ends), rules, Min)

1 ) 2 )

 12 Max-tree of p0's partial image. Max-tree of p1's partial image. Identify partitioned components. Utilize conflicts in canonical point labeling within the halo areas (area tuples). Tuples in this step are doubled, from large to small and inverted. Determine the edges of the boundary max-trees and represent them as tuples (root tuples).

  (a) Pan-sharpened, true-color image of Naples, Italy. (b) Central Milky Way, captured by the ESO, Chile.

Fig. 4 :

 4 Fig. 4: Benchmark images used in the experimental evaluation of the algorithm.

Fig. 5 :

 5 Fig.5: Execution time, speed-up and memory consumption curves of the proposed and Moschini's algorithm for increasing number of threads. The thread count for Moschini and Proposed (Shared) is increased locally on a single node. In case of the hybrid setting the number of threads and number of nodes are doubled alternately (refer to Table2).

Fig. 6 :

 6 Fig. 6: Execution time distribution of the proposed algorithm for the ESO 16bpp dataset. n signifies the number of utilized nodes and t the number of threads for on each node.

  Non-recursive version of Salembier's algorithm as presented by Carlinet et al.[START_REF] Carlinet | A Comparative Review of Component Tree Computation Algorithms[END_REF].

Wilkinson has proposed to use a priority queue and a stack, a combination of the algorithms of Salembier et al. and Hesselink et al., instead of using only a hierarchical queue for handling the pixel values during the flooding. Carlinet et al. [39] have proposed a non-recursive flooding algorithm variant of Salembier et al., which has strong Algorithm 1

while pqueue not empty do 34:

	15:	start pixel ← any point in Ω
	16:	PUSH(pqueue, start pixel)
	17:	PUSH(levroot, start pixel)
	18:	parent(start pixel) ← INQUEUE
	19:	⊲ 2. FLOODING:
	20:	loop
	21:	flood:
	22:	p ← TOP(pqueue); r ← TOP(levroot)
	23:	for all n ∈ N (p) such that parent(p) = -1 do
	24:	PUSH(pqueue, n)
	25:	parent(n) ← INQUEUE
	26:	if f (p) < f (n) then
	27:	PUSH(levroot, n)
	28:	goto flood	⊲ p is done
	29:	POP(pqueue)
	30:	parent(p) ← r
	31:	if p = r then INSERT FRONT(S, p)
	32:	⊲ 3. ROOT FIXING:
	33:

TABLE 1 :

 1 Overview of the worst-case time, space and message complexity of the distributed max-tree algorithm steps.

		Time	Space	Messages
	Image chunking Local max-tree	O( n p ) O( n p×t × log(d) + w × log(t))	O( n p ) O( n p + t × w)	O(1) -
	Root tuple generation O(w × d)	O(w × d)	-
	Area tuple generation O(w × log(p))	O(w)	O(log(p))
	Tuple resolution	O

TABLE 2 :

 2 Hybrid: multithreading+MPI.

	Nodes	1 1 2 2	4	4	8	8
	Threads 1 2 2 4	4	8	8	16
	Cores	1 2 4 8 16 32 64 128

TABLE 3 :

 3 Processing times (mean values in minutes and statistics) of the sequential Berger algorithm.

	Images	Mean µ StDev σ	CV	Min	Max
	Naples	13.54	0.484	0.001	13.53	13.55
	ESO 8bpp	113.76	6.011	0.001 113.62 113.91
	ESO 16bpp	184.56	9.877	0.001 184.35 184.73
	ESO 32bpp	185.67	19.457	0.002 185.34 186.21

TABLE 4 :

 4 Memory consumption (mean values in GB and statistics) for the different algorithms with Naples image when using a single node. For each threads setup, the average and the maximum resident set size of all the tasks in the job are reported, respectively.

		100.80	0.458	0.005 100.25 101.46	Average
		121.55	0.068	0.001 121.48 121.62 Maximum
	16	94.86 124.43	1.264 0.185	0.013 0.001 124.27 124.72 Maximum 93.58 96.30 Average
	24	94.27 124.82	1.714 0.241	0.018 0.002 124.46 125.05 Maximum 92.95 96.19 Average

Scientific and Grid Computing. His research interests include high productivity processing of big data in the context of scientific computing applications.

ACKNOWLEDGMENTS

The authors would like to thank Igancio Toledo and Martin Kornmesser for making the ESO/VVV Survey/D. Minniti image with the id eso1242a publicly available.


Liu Bo

Chen Guifen

The association rules algorithm based on clustering in mining research in corn yield

Keywords: corn yield, cluster analysis, association rules, frequent itemsets, Apriori Algorithm 0

With the popularization of agricultural information technology, the use of data mining techniques to analyze the impact of different types of soil nutrient content and yield of corn has become a hot topic in the field of agriculture. Association rule mining is an important part of the field in Data mining , association rules can be found associated with agricultural data attributes. This article will use cluster analysis and association rule to analysis correlation between corn yield and soil nutrient. Firstly compare different clustering algorithm to chooses the optimal algorithm, make data collected in scientific classification, and based on expert knowledge of the collected data into different levels; then determine the type and content of different soil by association rules corn yield and soil nutrient; final inspection algorithm is correct. The results showed that: comparing K-means, hierarchical clustering analysis, and PAM, K-means algorithm to determine the optimal clustering; K value can be determined at selected intervals.K is equal to3, 4 or 6,clustering effect is good according to Sil value when K from 3 to 10 . based on the principle of association rules, clustering algorithm to select a K value associated with the combination of rule 6; After clustering algorithm of association rules, support and credibility and improve degree of accuracy is better than not clustering; by mining association rules after clustering, a great influence on the different levels of soil nutrients in corn yield. The results for the corn yield provides intelligent decision support data.

Introduction

Corn is the world's most productive acreage and grain crops, it has fast growth rate and high yield characteristics [START_REF] Anthony | Plants and fertilisers as drivers of change in microbial community structure and function in soils[END_REF]. Moreover, with the development of animal husbandry and corn deep processing industry, corn has become the world's most important food crops, forage crops and cash crops [2]. 863 in support of national technology plan, the Shanghai, Beijing, Heilongjiang, Xinjiang, Jilin and other places to carry out tests to explore intelligent agriculture, the establishment of a number of precision agriculture experiment and demonstration area, and achieved gratifying results. China is a large agricultural country, China's grain output of corn yield very significant impact, where soil nutrients is one of the main factors affecting the yield of corn, so dig out the relationship between the different types and content of soil nutrients and corn yield is particularly important. But our data mining research later, there is no overall strength, compared with developed countries there is still a big gap, has seriously hampered the development of intelligent agriculture in China.

With the emergence and speed of data collection significantly improved, we urgently need new technologies and tools to the vast amounts of data into information and knowledge available to us. Data mining [START_REF] Hua | Different Fertilization Models oasis farmland soil microbial community abundance and Enzyme Activity[END_REF] is a large number of known data search and analysis to discover hidden potential in the data relationship, in order to predict the future. In the last decade, data mining as decision support has been rapid development [START_REF] Chu Now | Summary of data mining[END_REF]. Data mining has an important direction include cluster analysis [5] and association rule mining [START_REF] Xiang | A clustering algorithm of PSO & PAM[END_REF], the paper first by comparing the different clustering methods to select the optimal clustering algorithm, combined with association rules algorithm, dig out different relationship between soil nutrients and corn yield type and content. Data mining is used in agriculture to promote agricultural production and direction sustained, high-yield and effective means of security is important to protect the interests of farmers and national food security.

Clustering and association rules

Clustering is a way to simplify data through data modeling, which uses similarity to different data into different classes. So in the same class has a great similarity, in the different classes has a larger dissimilarity . Association rules are used to find links between things. Firstly, the concept and characteristics of each clustering algorithm is analyzed, based on the value of Sil select the optimal clustering algorithm; and then by the clustering algorithm is divided into different clusters to compare and choose the most appropriate number of clusters. Finally through the association rules to determine different kinds and content of soil nutrients and corn production relations.

Select Clustering Algorithm

Different algorithms have different characteristics and adaptability. For example, K-means algorithm [START_REF] Liu | Efficient Mining of Frequent Patterns Using Ascending Frequency Ordered Prefix-Tree[END_REF] When more dense clusters and cluster is obvious difference between clustering effect is good. But the noise and isolate sensitive data point. Specific steps are as follows: 1) Assign to each instance from its nearest cluster center to give K clusters; 2) were calculated for each cluster mean all instances, each of them as brand new cluster center. Repeat 1) and 2) until the position of the K cluster centers are fixed, assigned cluster is also fixed. Hierarchical clustering algorithm [START_REF] Aifen | Improved k-means initial clustering center selection method[END_REF] feature is not required prior to a given class number, the system can display the results in the clustering tree way, more suitable for data hierarchy, But to determine the distance matrix computation. Specific steps are as follows: 1) Each object is classified as a class, received a total of N classes, each class contains only one object. Distance between classes of objects they contain is the distance between. 2) find the closest two classes combined into one category, so the total number of class one less. 3) recalculate a new class clustering and all the old classes. 4) Repeat 2) and 3), until finally merged into a class so far.

PAM algorithm [START_REF] Li | Shuyun study courses relational structure [J] based on hierarchical clustering methods[END_REF] is less sensitive to noise and is not affected by the order of the input data, which is insufficient to determine the high computational cost of clustering centers required for large data clustering process is slow. Specific steps are as follows: 1) choose K objects as the initial cluster centers. 2) In addition to the open cluster center point of the sample to calculate the distance to each cluster center will classify the sample from the sample to the center of the nearest sample point. 3) and then calculated for each category, in addition to other classes outside the center of sample points and minimum distance to all other points, the minimum point as a new cluster centers. 4) Repeat 3) until the position of the two cluster centers unchanged.

Many data by understanding the structure and background information, you can know which algorithms are relatively good, but in many cases, our understanding of the data is not much, so we have to choose an objective criteria to evaluate [START_REF] Zhao | PAM clustering algorithm[END_REF]. Between the use of class compactness evaluated within the class separation and clustering is a common approach, one of the most classic is Silhouette indicators. Silhouette indicators both for the number of clusters optimal estimation can also be applied to evaluate the quality of clustering [START_REF] Harumi | Comparative data mining clustering technology research methods[END_REF]. Therefore, this paper Silhouette indicators as the number of selection and clustering clustering algorithm to determine the objective evaluation criteria.

Having provided a sample data set is divided into clusters, the clusters in the sample and the average of all other samples of dissimilarity or distance for the sample to an average of all samples of another class or dissimilarity distance, wherein, and. Thus, the index is calculated Silhouette samples are as follows:

(

All the samples in a cluster compactness of the average value of sil said tightness and separability; Sil average value of all samples of a data set may reflect the quality of clustering results , the greater the Sil value represents the better the quality of clustering.

The main steps of Silhouette clustering algorithm selection method based on the validity of indicators designed as follows: 1) setting a given data set, candidate cluster algorithm; 2) to specify the output number of classes K, respectively for each candidate cluster algorithm for data collection clustering; 3) calculating an average value of clustering results Sil each candidate clustering algorithm; 4) comparing the average value Sil, Sil value corresponding to the maximum average selected candidate clustering algorithm for the optimal algorithm.

In this paper, the national "863" plan "Maize Precise Operation System and Application" project demonstration base -Jilin Nong'an some experimental data from 2005 to 2010, and use MATLAB R2014a to cluster analysis. When K is equal to 3, the average Sil K-means algorithm is the highest; when K is equal to 4, the average value of K-means algorithm Sil also the highest. K-means algorithm proved optimal clustering algorithm, the results shown in Figure 1.

( ) ( ) ( ) [ ] ( ) ( ) { } t b t a t a t b t l , max Si - =

Clustering algorithm K value selection

Different values of K greater impact on clustering results, so traverse K 3-10 is the case, the corresponding value of Sil. as shown in picture 2. when K is equal to 3,4,6 clustering is better. According to the association rules algorithm is the core elements of frequent itemsets, therefore this article selects the K value is 6 k-means algorithm and its combined association rules algorithm.

Association Rules Algorithm

1993, R.Agrawal, who designed a program called Apriori algorithm [START_REF] Ying | A new algorithm for the optimal number of clusters[END_REF], the algorithm is the most influential algorithm, it laid the foundation for the Association Rules algorithm. This algorithm has two steps: The first step is mining frequent item sets. That support is greater than those specified by the user support selected projects, as frequent item set; the second step is based on frequent item sets to generate strong association rules. That association rules support and confidence are greater than or equal to the user specified support and confidence [START_REF] Cui | the Kingdom division An association rule mining for large transaction database algorithm[END_REF].

Support abbreviated as sup, Refers to a rule before or after a corresponding number of support as a percentage of total number of records. Formula is as follows:

.

(2)

Confidence abbreviated as conf, Said the DB contains the percentage of things at the same time also contains B. Formula is as follows:

.

(

Lift, it is the credibility of the ratio of the target density of the overall goal. So when the lift is greater than 1, the rules are better able to predict the results, rather than based on how often the data items to guess whether the result will be. Formula is as follows:

。 (4)

Case analysis

Based on the national "863" plan "corn precision operating systems research and application" project demonstration base -county agriculture of jilin province from 2005 to 2010, part of experimental data, experiment with clustering analysis and association rules method, validation and analysis of data mining method to the correlation between corn yield and soil nutrient.

Cluster analysis

In order to select frequent item sets, the first use of K-means algorithm on data collected scientific classification, data mining software Weka 3.7 data divided into six categories. And select the item appears set frequency greater than or equal to 15% based on expert recommendations [START_REF] Xiaoyu | An improved Apriori association rule mining algorithm[END_REF]. The third term due to the frequency appears set to 14%, the frequency of the sixth item appears set to 14%, it is discarded. The other four groups composed of frequent item sets. As shown in Figure 3.

( ) ( ) { } n B A TR TR B A B A È Ê = È = => sup sup ( ) ( ) ( ) A B A B sup sup A conf È = => ( ) ( ) ( ) ( ) ( ) B p A B p B B A conf B lift = => = => sup A

Data into Level

Based on the experience of experts of different data into different levels [START_REF] Xia | research on data mining[END_REF]. The N, P and K from A to F is divided into six different levels, the yield from A to C is divided into three different levels. And based on the knowledge of experts just four frequent itemsets obtained by clustering into different levels, as shown in FIG.

Relationship between the various soil nutrients and Yield

Respectively to observe the alkali solution N, P and K impact on production [START_REF] Zhaojun | under Different Fertilization film Maize dry matter accumulation and P uptake of[END_REF]. The results in FIG. 5 shown in FIG.

Data Association

To ignore an abnormal event and a small probability event. Selection criteria are: support between 15-100 percent, and lift the value of more than 1.0 and lift the top five values and yield-related association rules, as shown in Figure 6.

Examine the validity of

To test the accuracy of the algorithm, based on J48 algorithm [START_REF] Duga Silver | Save N Phosphorus stability control potassium fertilization on yield and nutrients accumulation Rice[END_REF] 10-fold cross-validation [START_REF] Yan | Chun decision tree algorithm for medical image data mining[END_REF]. The first data set is divided, which will in turn nine as training data, one as the test data, test. Each test will draw the appropriate accuracy. Then the average of 10 times as a result of the correct rate of accuracy of the method of estimation. As shown in Figure 7.

when yield is equal to A, the relationship between the various soil nutrients

A full-time to the Yield. Relationship between N, P and K between the three shown in Figure 8.

Conclusion and Analysis

By the above operation, we can see:

1. Choose a good clustering algorithm and K values of great influence on the clustering results.

Figure 1

 1 Figure 1 each algorithm of Sil value comparison

Figure 2

 2 Figure 2 different K value corresponding to the Sil value comparison

Figure 3

 3 Figure 3 clustering results and screening

Figure 4

 4 Figure 4 data into level

Figure 5

 5 Figure 5 is not clustering and clustering after the relationship between the soil nutrient and yield

Figure 6

 6 Figure 6 after clustering and clustering of relative contrast

Figure 7

 7 Figure 7 not clustering and clustering accuracy after contrast

Figure 8

 8 Figure 8 production for A while after not clustering and clustering of relations between the soil nutrient

Acknowledgment

Funds for this research was provided by National "863" project "corn precision operating systems research and demonstration" (2006AA10A309), National spark plan "corn" digital technology integration and demonstration(2008GA661003), Jilin provincial departments of world bank loan project(2011-Z20), Jilin province rural project in 2015(Study and demonstration of corn precision operation system based on Internet of things)

According to Comparative K-means, hierarchical clustering analysis, and PAM, K-means algorithm to determine the optimal clustering algorithm shown in Figure 1. 2. the value may be determined based on Sil K value chosen in the range from 3 to 10, when K is equal to 3,4,6 there will be a good clustering results, as shown in Fig. 2 3. After clustering accuracy of association rules algorithm is better than not clustering accuracy of association rules algorithm, show that using the association rules before clustering analysis of data processing is very necessary.. As shown in Figure 7. By the respective association rules diagram can be drawn: clustering support after the credibility and lift are better than not clustering, association rules show before using cluster analysis of data processing when the actual decision is significant . 4. the impact of different levels of soil nutrients for corn production is very high. When soil nutrients N content of C, P content of C, K content of E grade, grade A larger yield ratio.

Discussion and Outlook

In recent years, as we get the data and access to information ability increase [START_REF] Yingtang | Based on a fast cross-validation leave nuclear ELM online modeling[END_REF], expanding our database [START_REF] Hongwei | Monitoring data processing based on the state of the smart grid cloud[END_REF]. A growing number of agricultural production, scientific research, business management to use the database, and this trend will continue to rise. In this era of information explosion, information overload is everyone has to face the problem. How can you in a sea of information, find useful, hidden knowledge, as our top priority. Only make data resources, mining the data potential and useful knowledge, to make data have the effect of support decision making and prediction. Otherwise, a large amount of data it may become a waste, even burden [START_REF]Chen Hung Evaluation of Information Resources Quality[END_REF]. Therefore, data mining is a powerful guidance and direction to the development of the future, is a hot research field in the world today, its research has broad application prospects and great practical significance.