Mingyuan Jiu
email: mingyuan.jiu@liris.cnrs.fr

Christian Wolf
email: christian.wolf@liris.cnrs.fr

Graham Taylor
email: gwtaylor@uoguelph.ca

Atilla Baskurt
email: atilla.baskurt@liris.cnrs.fr

Human body part estimation from depth images via spatially-constrained deep learning

Keywords: Segmentation, Spatial layout, Deep learning, Convolutional networks, Depth images

Object recognition, human pose estimation and scene recognition are applications which are frequently solved through a decomposition into a collection of parts. The resulting local representation has significant advantages, especially in the case of occlusions and when the subject is non-rigid. Detection and recognition require modelling the appearance of the different object parts as well as their spatial layout. This representation has been particularly successful in body part estimation from depth images.

Integrating the spatial layout of parts may require the minimization of complex energy functions. This is prohibitive in most real world applications and therefore often omitted. However, ignoring the spatial layout puts all the burden on the classifier, whose only available information is local appearance. We propose a new method to integrate spatial layout into parts classification without costly pairwise terms during testing. Spatial relationships are exploited in the training algorithm, but not during testing. As with competing methods, the proposed method classifies pixels independently, which makes real-time processing possible. We show that training a classifier with spatial relationships increases generalization performance when compared to classical training minimizing classification error on the training set. We present an application to human body part estimation from depth images.

Introduction

Many computer vision problems can be solved in part by an initial step which segments an image, a video, or their constituent objects into regions, which are called parts in this context. The segmentation algorithm typically considers local appearance information, and frequently also models the spatial relationships between different parts. Unfortunately, considering these relationships within the segmentation process mostly amounts to solving constraint satisfaction problems or performing inference in a graphical model with cycles and a non sub-modular energy function, both of which are intractable in the general case. In this paper we address the problem of efficiently modeling spatial relationships without the need for solving complex combinatorial problems.

This general class of problems corresponds to various applications in computer vision. For example, pose estimation methods are also often naturally solved through a decomposition into body parts. A preliminary pixel classification step segments the object into body parts, from which joint positions can be estimated in a second step. The well-known system described in [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF], installed on millions of gaming consoles and taking as input images from consumer depth cameras, completely ignores spatial relationships between the object parts and puts all of the classification burden on the pixel-wise random forest classifier. To achieve its stateof-the-art level of performance, it required training on an extremely large training set of 2 • 10 9 examples.

A similar problem occurs in tasks where joint object recognition and segmentation is required. Layout

Z 3 i i N y i y 1 y 3 y 2 y 4 Z 2 Z 4 Z 1 Z i Z i y i (a) (b) (c)
Figure 1: Different ways to include spatial layout, or not, into learning part labels y i from features Z i for pixels i: (a) pixelwise independent classification, where spatial layout information is not taken into account; (b) A Markov random field with pairwise terms coding spatial constraints; (c) our method: pixelwise independent classification including spatial constraints N .

CRFs and their extensions model the object as a collection of local parts (patches or even individual pixels), which are related through an energy function [START_REF] Winn | The layout consistent random field for recognizing and segmenting partially occluded objects[END_REF]. However, unlike pictorial structures for object recognition [START_REF] Felzenszwalb | Pictorial Structures for Object Recognition[END_REF][START_REF] Felzenszwalb | Object detection with discriminatively trained part based models[END_REF], the energy function here contains cycles which makes minimization more complex, for instance through graph cuts techniques. Furthermore, the large number of labels makes the expansion move-type algorithms inefficient [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF]. In all cases, the underlying discrete optimization problem is very similar: an energy function encoding the spatial relationships in pairwise terms needs to be minimized. A typical dependency graph for this kind of problem is shown in Figure 1b: unary terms relate each label y i to a feature vector Z i , and pairwise terms encode prior knowledge on the possible configurations of neighboring labels y i and y j .

In this paper, we propose a method which segments an image or an object into parts through pixelwise classification, integrating the spatial layout of the part labels. Like methods which ignore the spatial layout, it is extremely fast as no additional step needs to be added to pixelwise classification and no energy minimization is necessary during testing. The (slight) additional computational load only concerns learning at an offline stage. The goal is not to compete with methods based on energy minimization, which is impossible through pixelwise classification only. Instead, we aim to improve the performance of pixelwise classification by using all of the available information during learning.

In each of the problems that we consider, the labels we aim to predict have spatial structure. Our proposed method uses an energy function to enforce a spatial con-sistency in learned features which reflects the spatial layout of labels. Unlike combinatorial methods, our energy function is minimized during training (i.e. while learning features) but is unused at test time. It is based on two main assumptions. First, different high-dimensional features with the the same label are embedded to a lower-dimensional manifold which preserves the original semantic meaning. Second is our belief that greater loss should be incurred when misclassification occurs between features coming from non-neighbor labels than features coming from the same or neighboring labels. In other words, the geometry of learned features, to some extent, reflects the spatial layout of labels. We will show that this new loss function increases the classification performance of the learned prediction model.

Another way of looking at our contribution is to interpret it as a way of structuring the prediction model of a learning machine. Classical techniques working on data represented in a vector space, like neural networks, SVMs, randomized decision trees, boosted classifiers, etc., are, in principle, capable of learning arbitrary complex decision functions if the underlying prediction model (architecture) is complex enough. In reality, the available amount of training data and computational resources available limit the complexity which can be learned. In most cases, only a limited amount of data is available with respect to the complexity of the problem. It is therefore often useful to impose some structure on the model. We already mentioned structured models based on energy minimization and their computational disadvantages. Manifold learning is another technique which assumes that the data, although represented in a high dimensional space, is distributed according to a lower dimensional manifold in that space. Semi-supervised learning uses a large amount of addi-tional training data, which is unlabeled, to help the learning machine better infer the structure of the decision function. In this work, we propose to use prior knowledge in the form of the spatial layout of the labels to add structure to the task of learning the decision function.

Another key aspect of our technique is end-to-end feature learning. The dominant methodology in computer vision, though changing in light of recent successes [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], is to extract engineered features such as SIFT [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF] or HOG [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF], pool responses, and learn a classifier from this fixed representation. Our objective is to apply learning at all stages of the pipeline, from pixels to labels. However, compared to contemporary Deep Learning approaches, we learn representations that are informed by the spatial structure of the part labels instead of simply their identity.

This paper proposes several contributions:

• We propose spatial learning for unsupervised pretraining of deep convolutional networks (i.e. learning all but the topmost layer). The features learned by spatial pre-training are more informative than classical features, as experiments with a non-spatial LR classifier show.

• We propose a framework which integrates spatial part layout into supervised learning of deep neural networks. We show that the gain of spatial learning is indeed obtained at extremely small cost: it improves the performance of the classifier with absolute zero increase in computational complexity of testing.

• To the best of our knowledge, we are the first to apply Deep Learning to the problem of body parts segmentation from depth images, obtaining promising results.

Related work

Our framework proposes to learn a feature extractor from raw data combining the spatial layout of labels, in order to produce a better decision function for segmentation. It is equivalent to learning a mapping function from high-dimensional space to a low-dimensional manifold space, therefore there is some partial overlap with dimensionality reduction. Unsupervised approaches for learning a mapping capturing global structure are well-known; most notably Principal Component Analysis (PCA) [START_REF] Jolliffe | Principal component analysis[END_REF] and Multi-Dimensional Scaling [START_REF] Cox | Multidimensional scaling[END_REF]. However, our aim is to embed based on the spatial layout of part labels, so we restrict our discussion to supervised methods. Neighborhood Components Analysis (NCA) [START_REF] Goldberger | Neighbourhood component analysis[END_REF] and its variants [START_REF] Salakhutdinov | Learning a nonlinear embedding by preserving class neighbourhood structure[END_REF] implicitly learn a mapping function from high dimensional space to low dimensional space while preserving the neighbourhood relationship defined by class labels. However, NCA is optimized for nearest neighbor classification and does not take into account structure within the labels, only their identity. Dr-LIM [START_REF] Hadsell | Dimensionality reduction by learning an invariant mapping[END_REF] is an online, non-probabilistic method which explicitly learns a non-linear invariant embedding function. Similarly, we parameterize our embedding with a convolutional neural network (ConvNet) [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF], however, like NCA, DrLIM does not consider structure in the labels. Our method differs from NCA and DrLIM by incorporating the spatial layout of labels into an energy function, rather than a binary notion of neighbors defined by datapoints with the same label.

ConvNets are a successful deep architecture widely studied in various applications, such as object recognition [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF][START_REF] Jarrett | What is the best multi-stage architecture for object recognition?[END_REF][START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], scene parsing [START_REF] Grangier | Deep convolutional networks for scene parsing[END_REF][START_REF] Farabet | Scene parsing with multiscale feature learning, purity trees, and optimal covers[END_REF] and connectomics [START_REF] Turaga | Convolutional networks can learn to generate affinity graphs for image segmentation[END_REF]. By parameter tying and feature pooling, ConvNets can automatically learn shift-invariant, discriminative low-and mid-level features from the raw pixels, avoiding the problem of generalization with handengineered features. [START_REF] Farabet | Scene parsing with multiscale feature learning, purity trees, and optimal covers[END_REF] have recently introduced multiscale ConvNets and applied them successfuly to Full Scene Labeling. Our work is similar in its motivation and the fact that we adopt their multi-scale approach to learn scale-invariant features. However, the way in which we approach the problem is very different. They apply sophisticated machinery based on optimal purity cover to search the best spatial grouping of feature vectors from which to predict a label. Our model has no such notion of adaptive pooling. Instead, we use weaklysupervised learning to introduce spatial context into the features. We believe the two approaches are complementary, and although beyond the scope of this paper, could be applied together.

Integrating spatial constraints into deep learning

In this section, we propose a method which learns to segment an image into parts which may be articulated. We pose this task as a typical classification problem, where each pixel of the image is classified into one of a discrete number of part labels L = {1 . . . L}. We aim to make decisions pixel-wise in order to minimize computational complexity.

We are given a set of M images {X1 , . . . , X M } and their associated labeled groundtruths. In our notation, the pixels of an image are indexed similarly index: X m = {X m (i)}. We seek to learn a segmentation model consisting of two parts:

• A parametric mapping function Z = f (X|θ f ) which
embeds each image X to a feature representation. This representation consists of Q maps, each of which have the same dimensions as X, and therefore can be indexed linearly:

Z(i) ∈ R Q .
The parameters θ f are learned from training data, taking into account Euclidean distances of pairs of features in the embedding space:

d(i, j) = ||Z(i) -Z(j)|| 2 (see Section 3.1).
• A classifier ŷ(i) = g(Z(i)|θ g ) which classifies the features Z(i) given trained parameters θ g giving an estimate ŷ(i) of the part label (see Section 3.2).

As is common in the Deep Learning literature, the embedding can be pre-trained in an unsupervised [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF] or supervised [START_REF] Salakhutdinov | Learning a nonlinear embedding by preserving class neighbourhood structure[END_REF] manner based on auto-encoding or some other kind of inductive principle. Then, the classifier and the embedding are jointly learned in a supervised way by minimizing some classification-based loss 1 . Our method proceeds in a similar way. We assume that the spatial part layout remains consistent across the different images of a corpus. In particular, adjacency information of parts is assumed not to vary across images. In body part estimation, for example, we suppose that the upper arm and the forearm are always adjacent. The next two subsections will describe how spatial constraints can be integrated into, respectively, the training procedure for the embedding f (•), as well as the training procedure for the classifier g(•). The two contributions can be applied independently, or combined. The result is a method proceeding as illustrated in Figure 1c: the information on the neighborhood layout is injected into a classifier working independently for each pixel i.

Spatial pre-training

Common ways of learning the embedding function f (•) are to minimize reconstruction error in an unsupervised auto-encoding setting [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF], or in a supervised setting, map pixels with the same labels to close-by points in feature space and pixels with different labels to distant points [START_REF] Salakhutdinov | Learning a nonlinear embedding by preserving class neighbourhood structure[END_REF][START_REF] Hadsell | Dimensionality reduction by learning an invariant mapping[END_REF][START_REF] Taylor | Learning invariance through imitation[END_REF].

Given the availability of ground-truth data (i.e. the part labels), we advocate its use in the pre-training step as supervised or semi-supervised training has been shown to produce a more discriminative embedding [START_REF] Salakhutdinov | Learning a nonlinear embedding by preserving class neighbourhood structure[END_REF]. However, instead of using classification loss as a training criteria, we aim to benefit from the knowledge of the neighborhood layout of the parts in the corpus. We introduce a new energy function which exploits the spatial layout of different parts:

E = i,j δ y(i),y(j) L S (i, j) + i,j ν y(i),y(j) L D (i, j) (1)
where y(i) and y(j) are the ground truth labels for pixels indexed by i and j, respectively; δ a,b is the Kronecker delta defined as δ a,b =1 if a=b and 0 otherwise; ν a,b is defined as ν a,b =1 if parts a and b are not neighbors in the corpus, and 0 otherwise. L S and L D are measured between pairs of features indexed by i and j in the feature space, which is embedded from raw data X through the embedding function. Note that we simply ignore contributions from pairs which have different labels but whose respective parts are neighbors in the corpus. As the exact boundaries between neighbors (given a human annotator) are somewhat arbitrary and often noisy, we choose to remain agnostic toward the relationship between features across a single boundary. We have chosen this scheme for simplicity, though we admit there are many other ways to design an energy function based on preserving spatial relationships. For example, if we could quantify the distance between labels (e.g. mean distance between body parts), such a "soft" neighbour criterion could be incorporated into the energy function. This is reserved for future work.

L S is the loss component for a pair of pixels having the same label, which has the effect of pushing their associated embeddings together. L D is a contrastive loss component for a pair having different labels which pulls their embeddings apart. Similar to [START_REF] Hadsell | Dimensionality reduction by learning an invariant mapping[END_REF], L S and L D are set to

L S (i, j) = 1 2 (d(i, j)) 2 , (2)
L D (i, j) = 1 2 [max(0, α -d(i, j)] 2 , (3)
where α > 0 is a margin. It controls the contrastive scope such that only dissimilar inputs which lie close-by in the embedded space contribute to the energy. The gradient of Equation ( 1) with respect to the parameters can be computed via the chain rule (backpropagation):

∂E ∂θ f = ∂E ∂Z(i) ∂Z(i) ∂θ f + ∂E ∂Z(j) ∂Z(j) ∂θ f (4)
where ∂E ∂Z(i) and ∂E ∂Z(j) can be obtained from ( 2) and (3).

∂Z(i)

∂θ f and ∂Z(j) ∂θ f depend on the particular architecture of the embedding function f (•) which is described in Section 3.3.

We finally want to point out that, although the energy function is defined on pairwise pixel terms, testing is done pixelwise, so testing does not require solving a combinatorial problem.

Supervised spatial LR learning

In Section 3.1, we introduced spatial constraints into unsupervised pre-training of deep learning architectures, which are ConvNets in our setting (see Section 3.3). In this section we show that a similar principle can also be applied to supervised learning of classifiers. In our work, we choose a ConvNet with a single fully-connected classification layer (i.e. the top layer is equivalent to logistic regression) which are trained end-to-end. The same principle, introducing spatial relationships into pixelwise classification, can also be applied to other classifiers, although the method may differ if the classifier is not learned to minimize classification error (or something akin to it, such as cross-entropy). In particular, in preliminary work we proposed a method to include spatial constraints into learning random forests by adapting the training algorithm based on maximizing gain in entropy [START_REF] Jiu | Integrating spatial layout of object parts into classification without pairwise terms -application to fast body parts estimation from depth images[END_REF].

Classification-based neural nets are typically trained to minimize cross-entropy. When the normalized outputs of the net are viewed as probabilities, this is equivalent to maximizing the log probability the net assigns to the true class. In the multi-class setting, this involves normalizing the outputs of the network via a softmax function and comparing them to the groundtruth label. However, minimizing cross-entropy does not take into account the layout of the part labels.

We propose the following new loss function, which is based on the ranking of class labels according to network output. For each input vector, a forward pass gives a network response for each class label, which can be used to rank the class labels in decreasing order. A loss can be defined based on the difference between this ranking and a desired target ranking, which is defined by the following properties:

• The highest ranked class label should be the target groundtruth label. This constraint is related to the entropy loss in traditional neural network learning;

• The next highest-ranked class labels should be neighbors of the groundtruth label in the class neighborhood definition of the corpus. We advocate that better generalization to unseen data can be achieved by forcing the net to learn these constraints.

An example for this is given in Figure 2, where the groundtruth label for the pixel is 1. The actual output ranks the groundtruth label at second place. The target ranking ranks groundtruth label 1 at first place, followed by labels 3 and 2 which, in this example, are neighbors of label 1.

Learning to rank is a classical problem in machine learning which has been addressed in the literature [START_REF] Burges | Learning to rank using gradient descent[END_REF][START_REF] Freund | An efficient boosting algorithm for combining preferences[END_REF][START_REF] Dekel | Log-linear models for label-ranking[END_REF]. We adopt a loss function similar in spirit to RankNet [START_REF] Burges | Learning to rank using gradient descent[END_REF], defined on pairwise constraints. Given a pair of labels (u, v), we denote by g u (Z(i)) and g v (Z(i)) the respective topmost (i.e. classification) layer outputs for location i and by o uv (i) = g u (Z(i)) -g v (Z(i)) their difference. The probability of label u being ranked higher than v is mapped through a logistic function:

P uv (i) =
e ouv(i) 1 + e ouv(i)

(5)

Given a target probability Puv (i), the loss function C uv (i) is the cross entropy loss [START_REF] Burges | Learning to rank using gradient descent[END_REF]:

C uv (i) = -Puv (i) log P uv (i)-1 -Puv (i) log (1 -P uv (i)) (6)
The target probability Puv (i) is set to λ>0.5 if class u is ranked higher in the desired ranking, and to 1-λ otherwise.

Given the properties of the desired ranking described above, the following two sets of pairwise constraints have been derived:

1. A set of L -1 constraints, where each constraint specifies that the groundtruth label is to be ranked higher than one of the other labels; 2. A set of constraints each one specifying that a label u, which is a neighbor of the groundtruth label, should be ranked higher than another label v, which is not a neighbor of the groundtruth label.

The loss function E for a single pixel is the sum over the pairwise constraints of the pairwise loss C uv :

E(i; l) = u =l C lu (i) + (u,v):ν u,l =0,ν v,l =1 C uv (i) (7)
where i is the index of the pixel, and we use the shortform l = y(i) as its ground truth label.

This loss function based on rankings provides a principled way of combining classification loss and spatial layout. Two types of constraints can be set through the λ parameter controlling the target probability Puv . Varying this parameter, different priorities could be given to different constraints.

It is important to note that this formulation allows one to include unlabeled data into the learning procedure in a semi-supervised setting. In this case, labeled images will be presented to the network with the loss function described above, whereas the loss function for unlabeled data does not contain the pairwise constraints including the groundtruth label -we consider neighbours of the network's strongest output instead of the ground truth neighbours.

As in Section 3.1, we would like to avoid confusion created by the term "pairwise", which here involves pairs of labels selected during training, which typically are of low cardinality compared to pixels. Testing proceeds independently for each pixel, conforming to the objectives of our work.

The deep learning architecture

While our mapping f (•) can be any differentiable architecture, we adopt a two-stage ConvNet [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF]. Each stage consists of a convolutional layer with a bank of filters, an element-wise non-linear squashing mapping (i.e. tanh function), followed by a spatial pooling and subsampling operator. A bank of feature maps are produced after convolving with the filters. The output of the first stage is fed as input to the convolutional layer at the second stage, where each feature map is connected to several maps in the previous stage. Connectivity is chosen uniform randomly before learning begins. The gradient of the energy function is computed with respect to the network output, and backpropagation is applied to update the parameters θ f . The parameters θ f consist of the filters and element-wise additive and multiplicative biases at the feature maps.

Many visual recognition problems require a large context due to complex interactions of different parts. However, the contextual size of ConvNets is limited by the filter size and sampling rate. Simply selecting larger filters does not address the fact that important cues may be observed at different scales. To overcome this dilemma, multiscale ConvNets [START_REF] Farabet | Scene parsing with multiscale feature learning, purity trees, and optimal covers[END_REF] are employed in our architecture, as shown in Figure 3. We use the notation f s (θ s f ), ∀s ∈ {1, . . . , N } to denote the output produced by each scale ConvNet where s indexes the scales and N is the total number of scales. By employing the weight sharing across scales, ConvNets can learn scale-invariant features.

A multiscale Laplacian pyramid X s , ∀s ∈ {1, . . . , N } is constructed for each image X in the pre-processing stage, where X 1 has the same size as the original image. Local contrast normalization is applied to the Laplacian images to ensure that the local neighborhood has zero mean and standard deviation. A batch consisting of pairs of patches2 is randomly extracted from each scale image X s . The patches are processed by the corresponding ConvNet f s , where the learned parameters θ s f are shared across the scales. The training procedure of the above architecture contains two steps. The first step is the spatial deep learning described in section 3.1, where the labels are only used to define spatial layout. We therefore call this stage weakly supervised pre-training. At this stage, the ConvNet parameters are initialized such that the features are consistent with the spatial arrangement of the labels. The second step is supervised spatial learning in section 3.2. A logistic regression layer parameterized with θ g is connected to the topmost feature maps of the ConvNet to predict the labels. We also apply a fine-tuning scheme in which the gradients from the LR are back-propagated to update the the ConvNet parameters θ f .

Experimental results

We experimented with our framework to validate spatial learning with multiscale ConvNets in different situations, with respect to the application of specific human body part estimation from depth images.

Unless otherwise specified, we use a 2-stage ConvNet in all of the experiments below. The first convolutional layer consists of 8 filters of size 9 × 9, followed by an element-wise tanh function and a 2 × 2 non-overlapping average pooling operator. The second convolutional layer consists of 32 filters of size 9 × 9, each of which combines 4 feature maps of the previous layer, with the same activation and pooling operator. The output is 32 feature maps. The receptive field for a single ConvNet on a single scale is 28 × 28. In a multiscale setting with a pyramid of 3 scales, with each scale consisting of 2 convolutional layers [START_REF] Farabet | Scene parsing with multiscale feature learning, purity trees, and optimal covers[END_REF], the receptive field is of size 112×112, thus a large spatial context is used to learn better features, comparable to best probe distance of 129 × 129 of the randomized forest in [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF].

Local contrast normalization is applied to the inputs of the ConvNet, and the pixels from the background or the context are set to an arbitrary high value (in our case, 4) to distinguish from the zero-distributed pixels of the object. The margin α for pre-training was set to 1.25 [START_REF] Hadsell | Dimensionality reduction by learning an invariant mapping[END_REF], and λ was set constant to 1 for all constraints in all experiments, i.e. Puv = 1 if u is ranked higher than v, otherwise 0. Weakly-supervised feature learning and supervised learning each used 30 epochs through the entire dataset, where mini-batch gradient descent is adopted. End-to-end fine-tuning, if applied, used another 30 epochs. Different learning hyperparameters (such as learning rates for different layers) were chosen empirically, e.g. 1 is set to 10 -6 for the first convolutional layer and 2 to 10 -5 for the second convolutional layer. When performing classification, we consider each pixel in isolation, applying its corresponding 32-dimensional feature vector as input to an LR. The parameters of the LR are shared at each pixel location. In the following experiments, we report mean pixel-wise accuracy.

The CDC4CV Poselets dataset (Holt et al., 2011) has been designed for human upper body part segmentation from depth images converted and rescaled to the grayscale domain (0-255). It contains 345 training and 347 test depth images of upper body poses shot with a Kinect module. All images have been pre-processed by background subtraction and annotated. The body parts are: head, neck, left shoulder, right shoulder, left upper arm, left forearm, right upper arm, right forearm, left hip, and right hip. Additionally, the part below the waist is given the label of type "other" (see the black area in Figure 5, column 2). We define a spatial layout of 28 neighborhood relations among the 11 parts calculated from a canonical pose where a subject stretches his or her arms, for example, head-neck, neck-left shoulder, Methods Accuracy Randomized forest [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF] 60.30% Spatial Randomized forest [START_REF] Jiu | Integrating spatial layout of object parts into classification without pairwise terms -application to fast body parts estimation from depth images[END_REF] 61.05% Single-scale (vanilla) ConvNet (LeCun et al., 1998) 47.17% Multi-scale ConvNet [START_REF] Farabet | Scene parsing with multiscale feature learning, purity trees, and optimal covers[END_REF] 62.54%

neck-right shoulder, left forearm-left upper arm, and so on. 900 pairs of patches were randomly chosen at three different scales for spatial pre-training.

For the application of part estimation from depth images, our baselines are the algorithms based on randomized decision forests described in [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF], as well as single-and multi-scale ConvNets with supervised end-to-end training without any spatial context. We also compare our framework with the recent spatial learning algorithm in randomized decision forests [START_REF] Jiu | Integrating spatial layout of object parts into classification without pairwise terms -application to fast body parts estimation from depth images[END_REF]. In our baseline experiments, a randomized forest of three trees with a depth of 9 is learned. For each node, the best offsets and threshold are chosen from 4000 offset candidates and 22 threshold candidates. More details can be found in [START_REF] Jiu | Integrating spatial layout of object parts into classification without pairwise terms -application to fast body parts estimation from depth images[END_REF]. From Table 1, it is clear that the large receptive field of the multiscale architecture is necessary to reliably estimate body part labels, and also that the Deep Learning-based architecture outperforms random forests on the CDC4CV dataset.

We investigate the performance of different combinations of our framework based on spatial pre-training of the purely convolutional layers (which we call ConvNet) and the classification layer (which we call LR). As a non-spatial baseline, we also implemented a "pixel-wise" version of DrLIM [START_REF] Hadsell | Dimensionality reduction by learning an invariant mapping[END_REF], a similar pretraining strategy in which spatial layout is not taken into consideration, i.e. δ a,b =1 if a=b and 0 otherwise; ν a,b =1 if a = b and 0 otherwise in equation 1. The results are shown in Table 2. For each setting (fine-tuning or no fine-tuning), spatial training outperforms non spatial training, and in many cases, the gains are high.

Examples of segmentation results are shown in Figure 5. From visual inspection we can see that the segmentation results produced by spatial learning are better than the ones by the non spatial methods. In particular, the segmentation produced by spatial learning is more consistent (less noisy), especially for the arms.

Spatial distribution of the error -The proposed method injects neighborhood relationships into the training algorithm. The question arises, whether the benefit of the proposed method only applies to pix- els at the boundaries, which would imply a rather trivial improvement. To rule this out, an in-depth analysis has been performed to verify that the proposed methods generally improve classification error independent of the pixel positions.

For each pixel in the test images the distance to the nearest part boundary was calculated using a distance transform, which allowed us to calculate histograms over these distances, where each bin corresponds to a range of distances. Histogram H g shows the number of pixels of the test images over distances, shown in Figure 4(a). This histogram largely depends on the sizes of the different body parts, and we can see that the pixel count decreases with distance. Distributions H s and H b are over pixels which are better classified by the baseline or by the proposed method, respectively:

• Distribution H s , shown in Figure 4(b), is the distribution of pixels which have been correctly classified by the proposed method but wrongly classified by the baseline method (DrLIM and classical LR with fine-tuning);

• Distribution H b , shown in Figure 4(c), is the distribution of pixels which have been correctly classified by the baseline method (DrLIM and classical LR with fine-tuning) but wrongly classified by the proposed method.

The normalized histogram H d , shown in Figure 4(d), illustrates the contribution of the proposed spatial learning method as a function of the pixels' distances from the part boundaries:

H d (i) = (H s (i) -H b (i))/H g (i). (8)
Our method gains over the very large majority of distances except for a few outliers. We see a strong bimodality, i.e. our method wins close-to and away from the borders. Moreover, the mode corresponding to pixels near the center of parts shows an even higher proportion of improvement than the mode closer to the border.

In particular, the classification rates are significantly improved between distance 50 and 60, which is more or less half of the receptive size. At the same time, this analysis verifies that our system indeed does benefit from the multi-scale framework.

Discussion

The experiments above demonstrate that our deep learning framework indeed benefits from the a priori additional information contained in the spatial part layout. Our framework can easily extend to semi-supervised learning, taking advantage of additional unlabeled training data. On this unlabeled data, the loss function would not include any terms based on classification loss (which require the ground truth part label), but only terms based on the spatial layout of the parts, as mentioned in Section 3.2.

It is worth noting that the two proposed contributions can be applied independently. In particular, the supervised spatial learning algorithm can be combined with any supervised or unsupervised training algorithms for ConvNets, for example auto-encoding or other algorithms reviewed in Section 2.

The quantitative evaluation measures presented in the different tables report pixelwise classification accuracy. This proves that the proposed learning algorithm improves pixelwise classification, although the injected additional information does not specifically state this. An additional advantage of our method is not captured by this evaluation: if a pixel is wrongly classified, then spatial learning will increase the probability that the wrong part label is a neighbor of the correct part label. In some applications this can be an advantage, for instance if the segmentation is post-processed to estimate joint positions, as in [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF].

Currently, we demonstrate one application of the proposed method to human part estimation in depth images. The approach is, however, not specific to depth images, and can be applied to other problems, as for instance object segmentation from RGB images, or semantic full scene labeling. In these cases, the context may be dynamic and the spatial layout might be more difficult to extract from the corpus, requiring learning techniques. This is reserved for future work.

The main computational complexity in our framework comes from the convolution operations, which can be significantly accelerated by parallel computation, e.g. on multi-core CPUs or GPUs. We implemented the system in MATLAB with multi-threaded convolutions provided by the GPU (CUDA) library (for training and testing), and the testing time by the IPP framework is also given as a comparison. We also give the running time of the randomized forest method implemented in C++. Although our training stage needs more computation time, testing per image is more efficient, since our method profits from a parallel implementation. Table 3 gives running times for testing and for training.

Conclusion

In this paper, we proposed a way to significantly improve classification performance in segmentation problems by integrating prior information on the spatial layout of image or object parts into a learning architecture. Compared to other spatial relationship learning algorithms, the energy function in our algorithm does not contain pairwise pixel terms, which makes it extremely fast. In particular, there is a zero increase in complexity compared to a standard pixelwise classifier of the same type in the test stage. We demonstrate applications of our approach to human body part labeling from depth images. Our algorithm outperforms the well known method based on randomized decision forests and other supervised feature mapping methods.

Figure 2 :

 2 Figure 2: The proposed loss function based on differences in ranking.

Figure 3 :

 3 Figure 3: Multiscale ConvNets framework[START_REF] Farabet | Scene parsing with multiscale feature learning, purity trees, and optimal covers[END_REF]. LCN means local contrast normalization. Each ConvNet contains several layers of convolutions and pooling described in 3.3.

Figure 4 :

 4 Figure 4: Histograms over different distances to part boundaries: (a) Hg -the histogram of all pixels; (b) Hs -the histogram of pixels for which the proposed method outperforms the baseline; (c) H b -the histogram of pixels for which the baseline outperforms the proposed method; (d) H d -is the normalized difference histogram.

Figure 5 :

 5 Figure 5: Classification examples from the CDC4CV dataset. (a) input depth image; (b) groundtruth segmentation; (c) appropriate baseline: randomized forest for CDC4CV; (d) DrLIM+LR without spatial learning; (e) our method (spatial pre-training and spatial LR learning).

Table 1 :

 1 Evaluation of different baselines on the CDC4CV dataset.

Table 2 :

 2 Results of different combinations of classical and spatial learning on the CDC4CV dataset. Fine-tuning means end-to-end training of the LR (top two layers) and ConvNet (remaining layers) with the same objective used to train the LR.

	Convolutional layers	LR	Fine-tuning Accuracy
	DrLIM	classical	no	35.10%
	DrLIM	spatial	no	41.05%
	spatial	classical	no	38.60%
	spatial	spatial	no	41.65%
	DrLIM	classical	yes	64.39%
	DrLIM	spatial	yes	65.12%
	spatial	classical	yes	65.18%
	spatial	spatial	yes	66.92%

Table 3 :

 3 Running times on the CDC4CV dataset of our proposed framework, compared with Randomized forest (RF). Training on 344 images includes pre-training, LR learning and fine-tuning , and testing time on one image is given.

	Task	Method Machine	Architecture	Time
	Testing	Ours	Xeon E56620, 8 cores, 2.4GHz MATLAB+IPP	00.41 sec
	Testing	Ours	Xeon E56620, 8 cores, 2.4GHz MATLAB+GPU	00.31 sec
	Testing	RF	Xeon E56620, 8 cores, 2.4GHz C++	00.73 sec
	Training Ours	Xeon E56620, 8 cores, 2.4GHz MATLAB+GPU 36.00 h
	Training RF	Xeon E56620, 8 cores, 2.4GHz C++	15.8 h

If pre-training is supervised, the common setting is to perform classification by k-nearest neighbor search in embedded space, requiring no additional learning step.

Since our loss function operates on pairs of features, i.e. particular locations i and j in the feature maps, we do not perform a full convolution while training. Instead, we extract patches corresponding to the receptive fields at locations i and j and map these patches to their respective feature vectors.


Florent Wachtel

Adrien Pilleboue

David Coeurjolly

Katherine Breeden

Gurprit Singh

Gaël Cathelin

Fernando De Goes

Mathieu Desbrun

Victor Ostromoukhov

Fast Tile-Based Adaptive Sampling with User-Specified Fourier Spectra

Keywords: Tiling Centroids Optimized centroids Sampling CR Categories: I.4.1 [Image Processing and Computer Vision]: Digitization and Image Capture-Sampling Sampling, Tile-Based Methods, Blue-Noise Distribution, General-Noise Distribution, Fourier Spectrum

Precomputation Runtime

: Using a carefully-designed set of subdivision rules, our approach constructs a self-similar, equi-area tiling that provides an even, unstructured distribution of tile centroids. Per-tile sampling points are then generated through offline optimization of the tile centroids in a series of local patches using any existing point set optimizer; the resulting offsets between original and optimized point locations are stored in a lookup table. These offset vectors are finally used at runtime to convert tile centroids into point set distributions with spectral properties closely matching those of the optimizer, but several orders of magnitude faster than current spectrum-controlled sampling methods.

Introduction

Sampling is key in a variety of graphics applications. In rendering, for instance, point sampling of continuous functions at strategically chosen locations allows for the accurate evaluation of spatial integrals. The spectral properties of these point samples have long been understood to be crucial in defining the notions of bias (systematic deviation from the ground truth), variance (signal-to-noise ratio), and aliasing (appearance of undesirable artifacts) in sampling-based evaluations. In particular, blue noise point sets [START_REF] Ulichney | Dithering with blue noise[END_REF]] are widely regarded as excellent due to their particular spectrum [START_REF] Pharr | Physically Based Rendering: From Theory to Implementation[END_REF].

Recent work has made more explicit the relationships between spectral content of a signal, spectral properties of sampling points, and resulting bias, variance, or aliasing [START_REF] Durand | A frequency analysis of Monte-Carlo and other numerical integration schemes[END_REF][START_REF] Subr | Fourier analysis of stochastic sampling strategies for assessing bias and variance in integration[END_REF]. Consequently, exercising precise control over the spectral properties of sampling points has been a recent subject of research [START_REF] Zhou | Point sampling with general noise spectrum[END_REF][START_REF] Öztireli | Analysis and synthesis of point distributions based on pair correlation[END_REF][START_REF] Gr Ünbaum | Blue noise sampling with controlled aliasing[END_REF]. Unfortunately, current algorithms offering spectral control are prohibitively slow-and thus poorly suited for what is arguably their most important application: Monte Carlo or Quasi-Monte Carlo integration. The only methods able to achieve the high throughput of blue noise sample points commonly needed in practical applications are based on precomputed tiles [START_REF] Lagae | Tile-based methods for interactive applications[END_REF]]. However, despite heavy offline optimizations, tile-based methods always exhibit spurious peaks and deficiencies in their associated Fourier spectra. These artifacts may be very pronounced if the tile barycenter is used as the sole sample per tile (Fig. 2), and remain visible even when there are multiple optimized samples per tile (Fig. 13) due to inherent patterns in the shape and positioning of the tiles.

Contributions. We present an efficient tile-based sampling system capable of generating high quality point distributions with controllable spectral properties at a rate of several millions of samples per second on commodity hardware-several order of magnitudes faster than current spectrum-controlled sampling methods. In sharp contrast to previous work, our approach relies on a deterministic set of production rules to produce self-similar, equi-area tiles whose centroids generate a blue noise distribution free of spurious spectral peaks. This tiling serves as a foundation for a lookup table based sampling system, in which the optimal positions of sampling points within each tile are determined by an offline optimization process to conform to a user-defined spectrum. Our approach thus allows for linear-time adaptive sampling of arbitrary density functions with precise spectral control over the resulting distributions.

To appear in TOG 33(4) SIGGRAPH 2014 irregular trihexes t (Sec. 3). The spectrum of irregular trihexes t is further improved via the optimization method described in Sec. 4.

Related work

In order to motivate and highlight our contributions, we briefly review the main approaches to generating point distributions with known desirable sampling properties and a few representative works that illustrate them best.

Stochastic sampling.

A first series of algorithms were formulated in the 1980s to eliminate aliasing artifacts in ray tracing [Dippé and Wold 1985;[START_REF] Cook | Stochastic sampling in computer graphics[END_REF][START_REF] Mitchell | Spectrally optimal sampling for distributed ray tracing[END_REF]. In particular, jittered sampling achieved a linear complexity compromise between uniform and fully-stochastic sampling. Poisson-disk distributions were also introduced via dart throwing [START_REF] Cook | Stochastic sampling in computer graphics[END_REF]], with many improvements and extensions proposed [START_REF] Mccool | Hierarchical Poisson disk sampling distributions[END_REF][START_REF] Ebeida | Efficient maximal Poisson-disk sampling[END_REF]. Current GPU implementations of stochastic sampling can generate over 100K samples per second [START_REF] Wei | Parallel Poisson disk sampling[END_REF][START_REF] Ebeida | Efficient maximal Poisson-disk sampling[END_REF]]. However, stochastic sampling methods such as jittering or Poisson-disk sampling produce distributions with relatively high variance by today's standards, and do not offer spectral control.

Optimal sampling. Over the last few years, higher quality sampling techniques have been proposed in which sample positions are spatially optimized to produce blue noise distributions [START_REF] Balzer | Capacityconstrained point distributions: A variant of Lloyd's method[END_REF][START_REF] Chen | Variational blue noise sampling[END_REF][START_REF] Schmaltz | Electrostatic halftoning[END_REF]Schlömer et al. 2011;[START_REF] Fattal | Blue-noise point sampling using kernel density model[END_REF][START_REF] De Goes | Blue noise through optimal transport[END_REF]. While these methods generally produce high-grade isotropic sampling, in practice their running times are prohibitive for large point sets. Recent methods [START_REF] Zhou | Point sampling with general noise spectrum[END_REF][START_REF] Öztireli | Analysis and synthesis of point distributions based on pair correlation[END_REF][START_REF] Gr Ünbaum | Blue noise sampling with controlled aliasing[END_REF]] have added control over the Fourier spectra of sample points by optimizing over a differential domain, an intermediate space between the spatial and Fourier domains [START_REF] Wei | Differential domain analysis for non-uniform sampling[END_REF]. Once again, these new methods suffer from prohibitive running times due to the need to iteratively proceed through inherently local optimizations.

Tile-based sampling. A decade ago, [START_REF] Cohen | Wang tiles for image and texture generation[END_REF] proposed a sampling method relying on the offline placement of samples within repeated tiles. Their approach was based on Wang tiles [1965] containing precomputed Poisson disk distributions, and was able to produce point sets containing desirable spectral properties with unmatched efficiency. Recursive tile subdivision [START_REF] Kopf | Recursive Wang tiles for real-time blue noise[END_REF]] and an improved management of tile borders [Lagae and Dutré 2006] furthered this work, resulting in fast and adaptive sampling of variable density functions. However, Wang-tile approaches all suffer serious drawbacks: the low count of prototiles and their placement on a square lattice induce a grid of peaks in the spectral domain, thus leading to sampling artifacts. A typical implementation of Wang tiles uses thousands of points per tile to alleviate these spectral peaks, at the cost of a significant increase in memory. Additionally, these methods rely on ranking the sampling points within each tile, bringing forth another critical flaw: for densities that are not a power of the subdivision factor, partial filling of prototiles is driven by a ranking algorithm [START_REF] Ulichney | The void-and-cluster method for dither array generation[END_REF]], not through offline optimization. This issue is clearly visible in the second row of Figure 11. Penrose tiles [START_REF] Ostromoukhov | Fast hierarchical importance sampling with blue noise properties[END_REF]] and polyominoes [START_REF] Ostromoukhov | Sampling with polyominoes[END_REF]] were proposed in lieu of Wang tiles as a convenient alternative.

With one sample per tile, these techniques allow for finer adaptivity when dealing with steep density functions. However, single-sample tiles result in rather poor Fourier spectra, as they exacerbate the presence of tiling structures which are overly simple or regular. While many other regular, semi-regular, or even irregular tilings may provide some improvement -likely at the price of implementation complexity -we are not aware of any self-similar tiling that fulfills the need for isotropic and even distributions.

Our Approach at a Glance

While existing tiling methods offer the best compromise between speed and quality of sampling, they fail to satisfy the two competing properties of using a deterministic, hierarchical construction, while generating a high quality isotropic, even, and non-periodic distribution of tiles. In this paper, such a novel tiling based sampling method is introduced. At its heart is an adaptive non-periodic tiling with equi-area tiles, paired with an offline computation of a table of optimized (spectrally-controlled) point sets defined over the tiles.

Rationale for our choice of tiling. From multiple failed attempts at improving previous tiling methods, we derived a few key insights that guided us in the design of our novel approach. First, we selected hexagons as the building blocks of our tiling rather than the squares used in Wang tiles [START_REF] Cohen | Wang tiles for image and texture generation[END_REF]] and polyominoes [START_REF] Ostromoukhov | Sampling with polyominoes[END_REF]; hexagons offer an additional symmetry axis, better intrinsic spectral properties, and only one adjacency relationship between neighboring blocks [Grünbaum and Shephard 1986].

Second, we derived a hierarchical and deterministic subdivision process based on trihexes (i.e., connected agglomerates of three hexagons, Fig. 2(e)) as prototiles to break the symmetries of the regular hexagonal lattice-much like polyominoes were unions of squares introduced to create more complex tiling structures on a square lattice. Finally, we devised a hierarchical, deterministic, and area-preserving border shuffling procedure that mimics a stochastic process in order to remove remaining frequency peaks in the resulting Fourier spectrum (Fig. 2(f)). We combine these three crucial features to create non-overlapping tiles whose barycenters form a high quality blue noise distribution. Note that each component is necessary: for instance, border shuffling without trihex subdivision only partially improves the spectrum of a tiling (Fig. 2(d)).

From tiles to controllable sampling.

Outline

We present our approach in four parts. Section 2 describes the creation of an irregular, equi-area tiling with hexagonal topology through border shuffling. Section 3 then introduces our trihex-based tiling, an extension of Ostromoukhov's polyominoes [2007], which, when used in conjunction with irregular hextiles, generates our tiling structure. The tiles' centroids will be shown to form a blue noise distribution free of spectral artifacts. We detail in Section 4 an adaptive sampling system with spectral control through an offline optimization of sample points on tiles. We finally demonstrate in Section 5 that our approach generates sampling points as fast as previous tile-based methods, but with state-of-the-art quality Fourier spectra on par with the best existing optimization methods.

Generating Irregular Hextiles

We begin with the construction of a hierarchical tilling in which each tile is an equi-area quasi-hexagon with randomized borders.

Quasi-Hexagonal Tiling

From the regular hexagonal lattice H (Fig. 3(a)) we create quasihexagonal regions we call hextiles, composed of λ hexagons from a finer regular hexagonal grid h (Fig. 3(b)): each hexagon Hi ∈ H is thus mapped to a hextile hi retaining the hexagonal structure of elements of H in that it still has six neighboring tiles. These hextiles form a tiling H of the plane and an equi-area partition of h. This construction defines a subdivision from H (spanned by the directions U1 and U2) to H (spanned by u1 and u2) whose subdivision factor is λ-chosen to be 37 to offer sufficient granularity (Appendix C). This subdivision can be recursively repeated down to any level. The unbounded, quasihexagonal lattice H over the finer regular hexagonal lattice h and generated by ( u1, u2). Each regular hexagon Hi in H corresponds, through subdivision, to a hextile hi made up of λ small hexagons.

Border Shuffling

Due to its regular structure, H and its subdivisions induce major spectral artifacts when used for sampling (Fig. 2(c)). To disrupt this structure, we introduce irregularity in the hextile borders while maintaining both the area and topology of the original hextiles {hi}, resulting in irregular hextiles { h i} that form, by construction, a tiling H of the plane.

✁ ✁ ✁

Discrete border alteration. We modify tile geometry using a stochastic process that reassigns hexagons of h along the original border from one hextile to an abutting one, thus maintaining the tiling property of the resulting irregular hextiles. To encode this border alteration, we use one bit per boundary hexagon to indicate if it is kept as is or reassigned to the adjacent hextile. One cannot pick random bits to alter shape: this would likely affect the evenness of the tiling by creating islands or holes, or by changing the tile area. Therefore, border bits must be carefully selected to guarantee the validity of the final irregular hextiling.

0 0 1 0 1 0 0 0 0 0 1 1 0 1 1 0 1 0 1 0 0 1 0 0 1 0 0 1 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 1 1 0 0 1 1 0 0 1 1 1 0 0 0 0 0 1 0 1 0 0 0 h i h k h j (b)
Figure 4: Effect of edge modifiers. From a set of hextiles hi and random valid edge modifiers (a), a set of irregular tiles h i in H is obtained (b) which maintains the topology of the initial tiling as well as each tile's initial area.

Irregular Hextiling

Given the choice of λ = 37, there is a total of 120 valid tripleedge modifiers, forming a sufficiently large set for our shuffling purposes. To construct an irregular tiling H, we begin with the regular tiling H and assign three independent triple-edge modifiers to each hextile as shown in Fig. 4(a). Reassigments are made based on the binary words bordering each regular tile hi, resulting in irregular tiles h i that are, by construction, also comprised of λ hexagons from h (Fig. 4(b)). This irregular hextile is guaranteed to be homeomorphic to a disk (Appendix A). Thanks to this areapreserving border shuffling, our subdivision process from H to H can thus be modified to instead map H to H, with the same subdivision factor λ.

Hierarchical Trihex Tiling

We adapt the polyomino tiling construction from [START_REF] Ostromoukhov | Sampling with polyominoes[END_REF]] to our hexagon-based approach, first without border shuffling. Unlike the original square-based method, we use prototiles that are m-polyhexes, i.e., connected conglomerates of equal numbers of hexagons [Grünbaum and Shephard 1986]. While an arbitrary value of m might be used, we choose m = 3 (i.e., trihexes) as it minimizes the memory requirement for our lookup tables while achieving our goal of improving spectral properties (Appendix C).

Prototiles

Subdivision Rules There are 11 trihex prototiles up to symmetries, shown in Fig. 5(a).

✁✂ ✄ ☎ ✄ ✆ ✝ ✝ t T (a) (b)
Moreover, we can extend the decomposition from H to H defined in Section 2.1 to trihexes: for each triplet of hexagons of H forming a prototile, we can trivially create a "quasi-trihex" on h, i.e., a group of three hextiles shaped like the original trihex (see Fig. 5(b)). In order to create a full hierarchical tiling, we decompose each quasi-trihex tile into a set of λ trihexes. This decomposition is performed offline using a simple discrete optimization as described in Appendix B.

The result is a cache of 11 subdivision rules {ti → Ti} transforming trihexes in H to λ trihexes on h, as illustrated in Fig. 5(b). This hierarchical tiling system is non-periodic, as follows from Statement 10.1.1 in [Grünbaum and Shephard 1986].

As substantiated by Fig. 2(e), these trihex subdivision rules {ti → Ti} form a tiling with a distribution of centroids whose Fourier spectrum has much reduced peaks when compared to the original spectrum of H. To further scramble this distribution, we now incorporate border shuffling.

Irregular Trihex-based Tiling

In order to mimic an ergodic stochastic process, border shuffling might be added at each level of the subdivision. However, this would defeat the purpose of using tiles for fast sampling, as it would create an exponentially large number of tiles. We thus must use our irregular tiling from Section 2 parsimoniously, adding just enough irregularity to the tiles without requiring undue computational and memory overhead. We achieve this balance by constructing a larger, but limited, number of subdivision rules.

Limited Irregularity via Bi-Periodicity We first limit the number of irregular hextiles in our tiling process through bi-periodicity in our triple-edge modifiers. We choose two generative integer vectors

g1 = 7 • u1-4 • u2 and g2 = 4 • u1+3 • u2
, where u1 and u2 are the vectors spanning our hextiling (see Fig. 3), such that the area A of the parallelogram formed by these two vectors is A( g1, g2) = λ • Ahex, where Ahex is the area of each hextile h 1 . The fundamental region of this bi-periodic construction is highlighted in Fig. 6 1 More generally, if λ is the n th centered hexagonal number, these vectors are chosen to be

g 1 = (2n+1) • u 1 -(n+1) • u 2 and g 2 = (n+1) • u 1 +n • u 2 . h i H g 1 g 2 u1 u2
Figure 6: Periodic tiling based on two integer basis vectors g1 and g2, creating a small set H of λ irregular hextiles { h i} (highlighted).

We then constrain our triple-edge modifiers to reflect the periodicity generated by these two vectors. Consequently, we must create a set of only λ different irregular hextiles { h i}i=1..λ := H repeated over H, as illustrated in Fig. 6. Notice that, by construction, we can naturally enumerate these irregular hextiles in H based on the (arbitrary) enumeration of the λ hexagons in h produced when subdividing a hexagon from H (see Fig. 3): one can therefore assign a hex-index ∈ [1..λ] per irregular hextile, uniquely determining the indices of the corresponding 6 adjacent hextiles. An example of this enumeration is given in Fig. 7.

Subdividing Irregular Hextiles. Each irregular hextile is further subdivided into smaller irregular hextiles from the limited set H: starting from an irregular tile h i with hex-indexed hexagons in h, we first apply one regular subdivision step. Then, if a hexagon h has hex-index i, the irregular tile used to replace h will be the i th irregular hextile in H. This subdivision step is illustrated in Fig. 7. The consistency of the subdivision system-in particular, the property that the irregular hextiles in Fig. 7 Ordering of Sub-Trihexes. Finally, trihexes in T are indexed by a rank (i.e., an integer in {1 . . . λ}). This choice of ordering is crucial to obtaining maximal quality during the adaptive sampling described in Section 4.3, as it assigns indices such that sequentially ranked trihexes are as physically distant from each other as possiblethus offering, between subdivision levels, a generalized Halton-like sequence of tile barycenters. We use the classical void-and-cluster method [START_REF] Ulichney | The void-and-cluster method for dither array generation[END_REF][START_REF] Kopf | Recursive Wang tiles for real-time blue noise[END_REF][START_REF] Ostromoukhov | Sampling with polyominoes[END_REF]] to find the ranking rank (t |t → T) of the irregular trihex t in each subdivision rule t → T. Rankings within the subdivision of a few trihexes can be seen in Fig. 10.

Concise Structural Indices

In order to efficiently access trihex lookup tables during sampling, we must be able to uniquely, yet concisely, identify a given tile at an arbitrary subdivision level. We address this important implementation issue by using hybrid structural indices as follows.

In our tiling T , an irregular trihex t k at subdivision level k could be uniquely defined by a sequence of subdivision rules from the library, i.e., as a sequence of trihex ranks2 :

rank (t 1| t 0 → T1) ⊕ . . . ⊕ rank (t k |t k-1 → T k ) , (1)
with t i a trihex in the set Ti. However, in order to support adaptive sampling (and thus, deep subdivision), k might become arbitrarily large, making trihex identifiers intractable in practice. We thus desire to formulate structural indices that are:

• finite, but from which trihexes of arbitrary resolution can be represented; • such that trihexes with the same structural index will share the same local configuration of trihexes, regardless of their subdivision level

These properties can be achieved by exploiting the cyclic structure in the sequence of trihex ranks, which leads to a limited number of possible local tile neighborhoods. In fact, we propose to use as structural indices a mix of neighborhood indices and ranks, thus offering a balance between the tile indexing proposed in [START_REF] Ostromoukhov | Fast hierarchical importance sampling with blue noise properties[END_REF]] (ranks only) and [START_REF] Ostromoukhov | Sampling with polyominoes[END_REF]] (neighborhood indices only).

Instead of explicitly representing the entire subdivision path from t 0 to t k , we characterize t k by the rank of its last l subdivision steps t k-l → T k , and replace the prepending path from t0 to t k-l with a neighborhood index based on the trihexes surrounding t k-l , which we denote n(t k-l ).

id(t k ) := n(t k-l ) ⊕ rank (t k |t k-l → T k ) . (2)
We found that using l = 1 provides a good balance between quality and lookup table size (see Appendix C). The neighborhood index n(t ) for a given irregular trihex t is a unique number representing the trihex adjacency around t , see Fig. 9. Such indices are computed from a local analysis of irregular trihex adjacencies in the subdivision rules: we observe that for each t , the number of unique neighboring trihex configurations n(t ) is 336 on average. An irregular trihex index is thus defined as one rank (a number between 1 and λ) and one neighboring index, which in our actual implementation leads to about 5M possible indices-thus 23 bits, however deep one needs to subdivide the tiles.

Adaptive Tile-based Sampling

At any subdivision level k, our tiling method generates an artifactfree blue noise distribution of trihex centroids with constant spatial sample density (Fig. 2(f )). However, single-sample-per-tile approaches produce high quality results only for powers-of-λ densities: in between two levels of subdivision, the spectral properties of the centroids' distribution are dictated solely by the ranking algorithm, resulting in suboptimal sampling. In order to remediate this issue and offer control over the distribution spectrum for the sampling of arbitrary constant densities, we now show how to generate spatiallyoptimized point samples for each tile in an offline stage, leading to high-quality adaptive sampling at runtime.

Lookup Sampling Table

We first construct a lookup table D to exert control over the spectra at intermediate levels of subdivision: for each tile, we will compute a set of optimized point samples at each range of ranks Rr := [0 . . . r] (with r < λ), where sample locations are stored as offsets from the tile centroid.

More precisely, we begin with a trihex for each neighborhood index and subdivide it several times, generating a large patch containing about N trihexes (the value of N is optimizer dependent; see Sec. 4.2). Inside each patch, each trihex t is subdivided once more, creating λ sub-trihexes. For each range of ranks Rr and for each trihex with structural index id in this patch, we create a point set containing the barycenter of the id trihex along with the N nearest sub-trihex barycenters with ranks in Rr (see Fig. 10). These points are then spatially optimized using any existing optimizer to improve their spectral distribution (discussed below). Finally, the resulting point set is converted into offsets from the barycenter of the id trihex, which are then stored in the lookup table entry D[r, id] for rank range Rr and structural index id. Note that our initial trihex subdivisions are such that we will run multiple optimizations of the same structural index. However, our construction of structural indices is such that two identical structural indices only share the same configuration of 1-ring of neighboring trihexes (see Fig. 9 in blue); farther neighboring trihexes may differ, but the offsets resulting from the optimizer are very similar in practice. Offsets from these multiple runs are thus averaged and stored in the lookup table as advocated in [START_REF] Ostromoukhov | Fast hierarchical importance sampling with blue noise properties[END_REF]].

(a) R 1 (b) R 2 (c) R 3 (d) R 37
Once the offsets for all structuring indices and all ranges of ranks are optimized separately, we convolve the results with a Gaussian kernel to smooth and enhance the offset coherence between rank ranges [START_REF] Ostromoukhov | Fast hierarchical importance sampling with blue noise properties[END_REF]]. The resulting lookup table, containing ∼ 90M elements (generated from ∼ 5M structural indices and λ = 37 ranges of ranks) is then able to handle adaptive sampling of arbitrary density functions as demonstrated in Fig. 11.

Optimizing Point Sets

Our approach can accommodate most existing point set optimizers that rely on local updates. Even methods based on global sample displacements such as [Schlömer et al. 2011] can be leveraged, at the cost of a postprocessing step to remap the final sample points to nearby trihex centroids. Based on the experiments described in Sec. 5.2, we found that our method produced high quality distributions even for very small patches around each possible irregular trihex (see the dotted circle in Fig. 9) when using a Lloyd based optimizer such as [START_REF] Balzer | Capacityconstrained point distributions: A variant of Lloyd's method[END_REF][START_REF] De Goes | Blue noise through optimal transport[END_REF]; we thus use a size N = 6 in this case. More general noise spectrum optimizers [START_REF] Zhou | Point sampling with general noise spectrum[END_REF][START_REF] Öztireli | Analysis and synthesis of point distributions based on pair correlation[END_REF][START_REF] Gr Ünbaum | Blue noise sampling with controlled aliasing[END_REF] warrant much larger patches, as they require accurate statistics on point distances. For these methods we use N ∼ 100, corresponding to a 5-ring vicinity. We are also able to tabulate our point sets to reproduce not only high-quality isotropic spectra, but also anisotropic spectra (see Fig. 14).

Adaptive Sampling

Upon completion of the offline optimizations and offset tabulation, we can now adaptively sample any importance function with the same approach described in previous tiling based methods [Ostro-moukhov et al. 2004;[START_REF] Ostromoukhov | Sampling with polyominoes[END_REF]]: from the importance function value, we locally select the most suitable subdivision level (through rounding) and a rank threshold to handle the density levels for octaves between power-of-λ levels. Samples are then generated by first performing an adaptive trihex-based irregular tiling (Section 3.2) at the given spatially-adapted subdivision levels, then placing, relative to the trihexes' centroids, precomputed samples from the lookup table based on the local rank threshold. The resulting samples (see, e.g., Fig. 12) exhibit the spectral distribution inherited from the point set optimizer that was used offline.

Results and Discussion

We now provide a series of tests and comparisons to demonstrate the high quality and interactive nature of our tile-based approach for various desirable spectra.

Testing against other Tile-based Methods

We begin by evaluating our results through comparisons with previous offline-optimized tile-based approaches-all of which are designed to produce blue noise point distributions.

Visual (spatial) quality. For constant sample density (Fig. 13), our approach stands out: sampling through Penrose tiles exhibits clear regularity and directionality, while polyominoes create distributions in which "holes" appear due to their relatively poor number of subdivision rules. Our results for density ramps (Fig. 11) also show fewer discrepancies in distances between neighbors through all density levels than previous approaches; in particular, Wang tiles lead to comparatively many holes and clusters.

Spectrum quality. Our tiling method is able to remove all peaks from the Fourier spectrum, even when using unoptimized tile centroids as the sampling points (Fig. 2(f)). After offline optimization, our one-sample-per-tile sampling has a blue noise spectrum that closely mimics the results of the best optimization techniques so far (Fig. 14(top)), and the anisotropy is both flat and at noise level.

In comparison, Fig. 13 shows that polyominoes retain visible residual peaks, while Wang tiles, even with 1024 samples per tile to improve their spectrum, exhibit anisotropy and a less flat low frequency region. For consistency, all power spectra, radial averages and anisotropies shown here were averaged over 10 distributions of 4K samples.

Speed. At runtime we can generate over 1M samples per second on commodity hardware, similar to most other tile-based methods.

As our approach is inherently parallelizable, one should be able to improve efficiency with specialized hardware.

Comparison with general-noise algorithms.

Spectrum shaping and Quality. Unlike previous tile-based sampling methods, we optimize our offsets offline in order to reproduce input-specified spectra with very high fidelity, as demonstrated in Fig. 14. The method of [START_REF] De Goes | Blue noise through optimal transport[END_REF] for blue noise and general-noise optimization techniques such as [START_REF] Zhou | Point sampling with general noise spectrum[END_REF][START_REF] Öztireli | Analysis and synthesis of point distributions based on pair correlation[END_REF][START_REF] Gr Ünbaum | Blue noise sampling with controlled aliasing[END_REF] were used successfully in our offline offset computations, and future point optimizers should be directly usable as well. We are also able to reproduce anisotropic spectra, even though our tiling is defined such that it forms an isotropic sampling: optimization of the offsets is sufficient to generate anisotropic distributions.

Speed. Again, we are able to generate millions of samples per second: our runtime complexity is independent of the target spectrum because our approach is based on lookup tables. Even the spectrum in Fig. 15, derived from observations of the spatial distribution of oak trees in a forest, can be adopted with no impact on timing.

Results of Adaptive Sampling

Our one-sample-per-tile output is able to capture varying density functions with smooth, blue-noise preserving transitions (Fig. 11).

Adaptive sampling is achieved by refining tiles to a depth proportional to the local density; our approach produces high quality results, as demonstrated in Fig. 12.

Generalization

While we focused on a particular implementation of our general approach that strikes a good balance between memory consumption and sampling quality, we note that different subdivision factors λ and/or other (m > 3)-polyhexes can be used without particular difficulty (see supplemental material). This flexibility can be leveraged to offer lower-memory or higher-quality alternative implementations when targeting more specific applications. For completeness, the supplemental material also provides an analysis of our stuctural index based on λ and on the number of ranks one may keep from the full indexing in Eq. 1.

Limitations

One drawback of our method is that it requires the selection of several parameters, which vary either the geometric complexity of tiles (λ, Sec. 2 and m, Sec. 3), or the extent of the irregular trihex neighborhood captured by the structural index (l, Sect. 3.3) or considered during the optimization step (N , Sect. 4.2). As discussed in Appendix C, these values represent a trade-off between point set quality and the size of the lookup table. We found that λ = 37, m = 3 and l = 1 provide the best compromise for this tiling system; but these values produce a lookup table containing 90M offset entries, potentially presenting a challenge for devices with limited memory. The last parameter, N (Sec. 4.2), does not affect the size of the lookup table, but does change the offline preprocessing time. This parameter varies from 6 to 100 depending on the optimization method the user selects to optimize point positions. Note that none of these parameters impact the timing of the final sample generation.

Another limitation is that this system is unsuitable for the generation of spatially structured distributions such as regular lattices, or of non-uniform distributions in the sense defined by [START_REF] Kuipers | Uniform Distribution of Sequences[END_REF]. In other words, our method is best suited for distributions that are locally characterized, and thus generate bounded offsets. Consequently, we follow the work of Heck and colleagues [2013] and focus on spectra from the blue noise family.

Conclusion

Our sampling approach provides an extremely effective way to generate high-quality spectrum-specified sample distributions. From a novel tiling with blue noise properties, precomputed offsets with user-defined spectral properties are used to sample any density function in linear complexity at a rate of millions of samples per second on commodity hardware. While our new tiling procedure exhibits sampling properties that are far superior to current methods, we believe our approach can be further improved upon. First, our lookup tables have a non-negligible memory footprint which may limit implementations on GPUs. Furthermore, it would be desirable to have multiple spectral targets embedded in our hierarchical tile-based sampling system. The ideas presented in [START_REF] Belcour | 5D covariance tracing for efficient defocus and motion blur[END_REF]] may be pertinent to this extension.

Acknowledgements. This project was supported by the ANR excellence chair (ANR-10-CEXC-002-01) and digitalSnow program (ANR-11-BS02-009), as well as NSF grants CCF-1111943 and CCF-1011944, and a Google graduate fellowship. The authors are grateful to the anonymous reviewers for their constructive comments and suggestions to improve the final version of this paper. We also thank Kartic Subr for helpful discussions, and Patrick Mullen for his insightful comments.

A Valid triple-edge modifiers

The two sufficient conditions on the validity of a triple-edge modifier explained in Section 2.2 are trivial to prove. The first requirementhaving the same number of non-zero bits in the modifier words A,B, and C-enforces that each tile preserves its original area. Indeed, following Fig. 4(a), any 1 in A means that a hexagon of hj will be added to hi, which implies that there must also be a corresponding 1 in B and C, indicating that a hexagon of h k will be added to hj and another moved from hi to h k . Thusly, all modified hextiles will still contain λ hexagons. The second requirement-having none of the words end with "01"-is easily shown to imply homeomorphic shuffling (i.e., no islands or holes) through enumeration of all configurations of A with at least one non-zero entry (dark orange hexagons):

Grey hexagons are those in hi that may be removed from hi depending on surrounding triple-edge modifiers. We see that in all but the rightmost case (word ending with "01"), every added hexagon is adjacent to a hexagon of hi that will remain in h i (light orange hexagons), thus resulting in a modified tile h i homemorphic to hi. In our supplementary material, we provide a formal proof that these two conditions are sufficient for arbitrary centered hexagonal numbers.

B Decomposition of T into irregular trihexes

As described in Sec. 3.2, an important step of our approach consists of decomposing T into irregular trihexes from T . While this is likely to be an NP-complete problem, we purposely kept small both the cardinality | T | = 11 • λ and the number λ of irregular trihexes t i within each element. We can thus directly explore the combinatorial search space to get a solution. In fact, one can design a very efficient algorithm to quickly construct a wide set of solutions by using the concept of discrete capacity-constrained Voronoi tessellation (CCVT) proposed in [START_REF] Balzer | Capacityconstrained point distributions: A variant of Lloyd's method[END_REF]]. Indeed, we aim to decompose T into cells (defined as sets of irregular tiles h ) with capacity 3 (each h having a capacity of one), and we know that the tessellation induces a valid decomposition of T if each cell defines an irregular trihex (i.e., if its three irregular tiles are adjacent). Following [START_REF] Balzer | Capacityconstrained point distributions: A variant of Lloyd's method[END_REF]], we first attach random labels {1 . . . λ} to each tile and perform label swaps if the proposed change minimizes the CCVT energy. After several swaps, we check if the resulting decomposition is valid-that is, if each cell is an irregular trihex. If not, we perform local random label swaps and continue iterating. This minimization of the CCVT energy allows us to quickly and robustly obtain valid and non trivial decompositions in only a few iterations. Note that this process is so fast to converge that one can generate multiple valid decompositions of T and pick the decomposition for which a maximum number of different prototiles is used: this process further improves the mixing of prototiles within the decomposition.

C Parameter Selection

Here we provide a brief review of the parameters used in the results shown. For a more complete treatment, along with additional illustrations and spectra, we refer the reader to the supplementary Polyhex size, m: Smaller polyhexes (m = 1, see Fig. 2, or m = 2) produce inferior blue-noise spectra. Larger polyhexes lead to more possible tile configurations, but also more neighborhood indices. The greatest jump in quality is found when increasing m from 2 to 3, with the advantages for larger polyhexes outweighed by the corresponding increases in lookup table size.

Structural index depth, l: Increasing this parameter expands the neighborhood referred to by each structural index, with an expo-nential effect on the number of structural indices generated. Vast improvements to sample distribution and spectral quality are observed when increasing l from 0 to 1; the corresponding effects when increasing from 1 to 2 are limited. Thus, l = 1 is found to adequately characterize tile neighborhoods at an acceptable cost. our tiling based method using each of these optimization methods for our offline offset computations. Each example includes Fourier spectrum, point distribution, radial average and anisotropy of power spectrum. Notice the very limited degradation of the spectrum generated by our tiling approach (see also Fig. 13).

Figure 2 :

 2 Figure 2: Spectral characteristics of various tile-based sampling methods for one sample point per tile. Top row: tiling and tile centroids. Bottom row: Fourier spectra of tile centroids. Circles on the right portion of each spectrum indicate peaks of size 0.1 and greater (DC value of 1). Left to right: (a) Penrose tiles, (b) polyominoes, (c) regular hexagons, (d) irregular hextiles h (Sec. 2.2), (e) regular trihexes t (Sec. 3.1), (f)

Figure 3 :

 3 Figure 3: Quasi-hexagonal tiles. (a) The unbounded, regular hexagonal lattice H generated by ( U1, U2). (b) The unbounded, quasihexagonal lattice H over the finer regular hexagonal lattice h and generated by ( u1, u2). Each regular hexagon Hi in H corresponds,

  modifiers. A simple characterization of valid border bits is achieved by considering triplets of connected border tiles. For each, we associate a triplet of binary words A := a0a1a2a3, B := b0b1b2b3, and C := c0c1c2c3 called a triple-edge modifier using the convention depicted in the inset figure, where 0 indicates no local reassignment. As discussed in Appendix A, each triple-edge modifier must consist of three words A, B, and C containing an equal number of nonzero entries, none of which terminates in the sequence "01." Note that this straightforward characterization of valid triple-edge modifiers is in fact sufficient for any λ.

Figure 5 :

 5 Figure 5: (a) 11 configurations of trihexes on H. (b) The library of subdivision rules {t → T}, where the quasi-trihex tile associated with t after subdivision is partitioned into a set of λ = 37 trihexes.

Figure 7 :Figure 8 :

 78 Figure7: From an irregular hextile h i with indexed hexagons (left), we apply one regular subdivision step to construct sub-hextiles hi (center). Replacing each of these λ sub-hextiles with a different irregular hextile from H creates a seamlessly subdivided irregular hextile (right).

Figure 9 :

 9 Figure 9: Neighborhood of a trihex. Two trihexes (dark blue) with identical neighborhood indices have the same 1-ring of neighboring trihexes (light blue) but may have different (n > 1)-rings. Point-set optimizations start with the centroid of each tile and are averaged over all possible configurations of a trihex's vicinity (dotted circle).

Figure 10 :

 10 Figure 10: Optimization snapshots for various sample densities (i.e., various rank ranges Rr). Black dots mark the centroids of sub-tiles (colored); arrows indicate the sampling offsets computed in the first iteration of the optimization loop.

Figure 11 :

 11 Figure 11: Adaptive sampling of a quadratic density ramp with 2000 points for various sampling algorithms.

Figure 12 :

 12 Figure 12: Importance sampling for rendering light integration; high-dynamic range image courtesy of Paul Debevec.

FrequencyFigure 13 :

 13 Figure 13: Comparison with other tile-based sampling methods for constant density function. Top to bottom: power spectra, sample distribution, radially-averaged power spectra and radial anisotropy. Notice the dramatic improvement in blue noise profile and isotropy.

FrequencyFigure 14 :

 14 Figure14: Optimization-based vs. tile-based sampling methods. Left: results of blue noise and general noise optimization methods. Right: our tiling based method using each of these optimization methods for our offline offset computations. Each example includes Fourier spectrum, point distribution, radial average and anisotropy of power spectrum. Notice the very limited degradation of the spectrum generated by our tiling approach (see also Fig.13).

FrequencyFigure 15 :

 15 Figure 15: Distribution (left) of oak trees based on a measured distribution from [Pommerening 2002]. Visualization (center and right) of the resulting forest.

Table 1 :

 1 Throughout this paper, we will denote by H the regular hexagonal lattice and h a scaled down version of this lattice. Hexagons Hi ∈ H are subdivided into quasi-hexagonal tiles hi we call hextiles, formed by the union of λ smaller hexagons from h. Here, λ is a centered hexagonal number Cn, with Cn:= Main notations used in the description of our approach.

		As hinted by the subtle
	anisotropic structure visible in Fig 2(f), the removal of spectral peaks
	through our novel tiling may not be sufficient for most sampling
	purposes. Equipped with a non-periodic and even distribution, we
	describe the construction of lookup tables that associate each tile
	with an optimized point location, removing residual anisotropy and
	allowing for the efficient generation of high quality samplings at
	runtime. Existing optimization procedures may be used in this
	offline step to generate samplings with given spectral properties.
	Point sets are stored as offsets from the centroid of each polyhex
	and are made coherent across subdivision levels, thus guaranteeing
	high quality adaptive sampling of arbitrary density functions.
	Notation. 1, 7, 19, 37, 61, ... for n = 1, 2, 3, 4, ... [Conway and Guy 1996];
	while we adopt λ := C4 = 37 to minimize memory requirements
	and complexity, other Cn can be used as subdivision factors (Ap-
	pendix C). The partition formed by these hextiles is denoted as H.
	Hextiles are further made irregular by randomly reassigning small hexagons on h to adjacent hextiles while preserving each hextile's
	connectedness and area, creating a tiling H of irregular hextiles { h i}. We use trihexes (i.e., connected agglomerates ti of three
	hexagons) as prototiles, such that irregular trihexes (i.e., connected
	agglomerates t i of three irregular hextiles) will form the 407 tiles,
	unique up to rotations and translations, of our final tiling T (Sect.
	3.2). Table 1 illustrates our notations used throughout this paper.
	Notation	Description	See also
	Hi Regular hexagon from regular grid H. Sec. 2.1, Fig. 3
	hi Hextile; partition of h (smaller regular grid) into hextiles is denoted H.	Sec. 2.1, Fig. 3
	h i	Irregular hextile; partition of h into irregular hextiles is denoted H.	Sec. 2.2, Fig. 4
	ti Regular trihex on H.	Sec. 3.1, Fig. 5
	t i	Irregular trihex; partition of h into ir-regular trihexes forms tiling T .	Sec. 3.2, Fig. 8

Table 2 :

 2 Top to bottom: power spectra, sample distribution, radially-averaged power spectra and radial anisotropy. Notice the dramatic improvement in blue noise profile and isotropy. Overview of the properties and timings of sampling algorithms. Timing were measured on a single core 3.2 GHz CPU, except those marked with an asterisk *, which utilized an NVIDIA GTX 660 GPU. materials available on the project website.Subdivision factor, λ: This parameter determines the number of hexagons into which tiles are subdivided, and is thus constrained to the centered hexagonal numbers (1, 7, 19, 37, 61, ...). Smaller λ yield tile sets with insufficient variety to overcome the limitations common to previous tile-based sampling methods (Sec. 1). On the other hand, larger values lead to finer subdivision, reducing the influence of each shuffled border hexagon (Sec. 2.2) on its tile centroid and increasing the lookup table size. Using λ = 37 limits lookup table size, provides sufficient tile variety, and maximizes effectiveness of border shuffling.

			Qualitative Analysis			Timing (seconds)
	Method	Speed	Memory footprint	Distribution Quality	Spectral Control	10k	50k	1M	10M
	[Schlömer et al. 2011]	very poor	good	good		157.2	4880	-	-
	[de Goes et al. 2012]	poor	good	excellent		6.135	65.00	10738	-
	General Noise	poor	good	good		10.11* 13.91* 147.5*	-
	[Ostromoukhov et al. 2004]	excellent	satisfactory satisfactory		0.730	4.55	11.8	80.14
	[Kopf et al. 2006]	excellent	poor	good		0.004	0.018	0.345	3.42
	[Ostromoukhov 2007]	excellent	poor	good		0.053	0.254	18.9	19.6
	Ours (blue noise target)	excellent	poor	excellent		0.016	0.067	0.989	9.28
	Ours (general noise target)	excellent	poor	good		0.016	0.067	0.989	9.28

Trihex-based Tiling via Irregular HextilesOnce the border shuffling procedure is applied, the resulting tiling H will contain a rich set of irregular hextiles. While improved, its underlying quasi-hexagonal structure still results in an unsatisfactory distribution of tile centroids (Fig.2(d)). This issue is addressed next through the use of m-polyhex tiles (here with m = 3, see Appendix C), extending the polyominoes method[START_REF] Ostromoukhov | Sampling with polyominoes[END_REF]] to our irregular hextiling setup. The resulting construction of subdivision rules based on irregular trihexes will serve as the basis of the tiling T used in our sampling approach.

If ranks are integer numbers on b bits, ⊕ is the concatenation operator on bits. The trihex identifier would thus require k • b bits.