Junlin Ouyang

Gouenou Coatrieux

Ph. D Huazhong Shu

Robust hashing for image authentication using quaternion discrete Fourier transform and log-polar transform

Keywords:

Robust hashing for image authentication using quaternion discrete Fourier transform and log-polar transform

Junlin Ouyang, Gouenou Coatrieux, Huazhong

INTRODUCTION

With the widespread use of sophisticated image editing tools, image contents can be easily tampered or forged.

Verifying image authenticity is therefore a major issue in many applications. Robust image hashing is widely applied in image authentication [START_REF] Ahmed | A secure and robust hash-based scheme for image authentication[END_REF][START_REF] Lei | Robust image hash in Radon transform domain for authentication[END_REF][START_REF] Zhao | Robust hashing for image authentication using Zernike moments and local features[END_REF]. Schematically, image hashing methods extract essential image features from which a short binary or real number sequence, called hash, is generated to represent the image content. Such a hash should be robust to image content-preserving operations while being sensitive to malicious tampering ones.

Because robust image hashing captures the main image characteristics, it has attracted interest for other applications like image forensic [START_REF] Lu | Multimedia forensic hash based on visual words[END_REF][START_REF] Battiato | A robust forensic hash component for image alignment[END_REF], image retrieval [START_REF] Liu | Spline regression hashing for fast image search[END_REF][START_REF] Zhu | Sparse hashing for fast multimedia search[END_REF], digital watermarking [START_REF] Fridrich | Robust hash functions for digital watermarking[END_REF][START_REF] Cannons | Design and statistical analysis of a hash-aided image watermarking system[END_REF], and so on. Regarding image authentication, a robust image hash should also have good anti-collision (discriminative) capability for visually distinct images as well as a satisfactory level of security in order to make very difficult for an adversary to forge the hash value. To meet these requirements simultaneously, the construction of a robust image hashing is still a challenging task.

In general, the construction of an image hash is based on three basic steps, i.e., pre-processing, image feature extraction and construction of hash. Among them, the most critical one is certainly the image feature extraction step [START_REF] Swaminathan | Robust and secure image hashing[END_REF]. Existing feature extraction methods for robust image hashing can be roughly classified into the following categories. [START_REF] Fridrich | Robust hash functions for digital watermarking[END_REF] used the DCT to capture the essential features of image blocks. They observed that it is very difficult to change the correlation of the low-frequency DCT coefficients without tampering the content of an image. Therefore, low-frequency DCT coefficients can be utilized as features to build the hash. De Roover et al. [START_REF] Roover | Robust video hashing based on radial projections of key frames[END_REF] defined a radial variance vector (RAV) based on the radial projection of the image pixels, and then applied the DCT to compress the RAV feature vector and construct the image hash. This method is robust to content-preserving operations and small angle rotations.

Discrete Cosine Transform (DCT)-based image hashing methods: Fridrich and Goljan

Discrete wavelet transform (DWT) -based image hashing methods: Ahmed et al. [START_REF] Ahmed | A secure and robust hash-based scheme for image authentication[END_REF] used a wavelet transform to extract the image features. Since the wavelet transform has a good time-frequency localization property, their method can locate tampered regions with a good accuracy but at the price of a longer image hash. Wu et al. [START_REF] Wu | A novel image hash algorithm resistant to print-scan[END_REF] combined the Radon transform (RT) and the DWT to deal with print-scan attacks. Tang et al. [START_REF] Tang | Perceptual image hashing using local entropies and DWT[END_REF] developed an image hashing scheme after observing that the entropy of pixel blocks, a measure used to characterize image texture, is approximately linearly changed after content-preserving operations. Then, they applied DWT to image block's entropies to realize feature compression and construct the image hash. Recently, color vector angle combined with DWT [START_REF] Tang | Robust image hashing via colour vector angles and discrete wavelet transform[END_REF] were applied to robust image hashing. Their results show good robustness to common content-preserving operations and small angle rotation. Liu et al. [START_REF] Liu | Wave atom transform generated strong image hashing scheme[END_REF] utilized the wave atom transform to extract image features arguing that this approach has a sparse expansion and is capable to better capture texture properties. Furthermore, they observed that the coefficients of the third scale band are more suitable to serve as image hash features than the other scale bands.

Discrete Fourier transforms (DFT)-based image hashing methods:

Most DFT based approaches are combined with other transforms in order to resist to geometric distortions. For example, Swaminathan et al. [START_REF] Swaminathan | Robust and secure image hashing[END_REF] proposed a robust image hashing scheme based on the Fourier-Mellin transform. Qin et al. [START_REF] Qin | Robust image hashing using non-uniform sampling in discrete Fourier domain[END_REF] introduced another scheme where a secondary image is first obtained after a rotation projection similar to the RAV of the image. A non-uniform sampling is then performed to extract the image features after applying the DFT. This method is robust to small angle rotations. Lei et al. [START_REF] Roover | Robust video hashing based on radial projections of key frames[END_REF] first carried out a RT and then computed moment features before applying DFT on these moments. The first fifteen significant DFT coefficients were then normalized and quantized to obtain the image hash value. This method shows satisfactory results when facing geometrical distortions.

Matrix decomposition-based image hashing methods: Kozat et al. [START_REF] Kozat | Robust perceptual image hashing via matrix invariants[END_REF] used singular value decomposition (SVD) to get robust image features and to generate an image hash. Their hashing algorithm consists of two major steps. In the first one, intermediate features are extracted from pseudo-random (PR) semi-global regions via SVD. In the second step, the SVD is again applied to the intermediate features so as to construct the final hash. Monga et al. [START_REF] Monga | Robust and secure image hashing via non-negative matrix factorizations[END_REF] also introduced a new PR signal representation method using non-negative matrix factorization (NMF) to capture the features and to form the image hash. Their method shows better performance than the SVD based method.

Recently, Tang et al. [START_REF] Tang | Robust perceptual image hashing based on ring partition and NMF[END_REF] proposed a robust image hashing method based on ring partition and NMF. Their results show good performance to common content-preserving operations and large angle rotation operation.

Others image hashing methods:

In [START_REF] Xiang | Histogram-based image hashing scheme robust against geometric deformations[END_REF], Xiang et al. proposed a histogram-based image hashing scheme, which is robust to geometric distortions but not to additive noise, brightness adjustment and contrast enhancement. Battiato et al. [START_REF] Battiato | Robust image alignment for tampering detection[END_REF] adopted an image representation based on a set of SIFT features (called "bag of features", BOF) to construct the hash and explored a non-uniform quantization of histograms of oriented gradients (HOG) to get tamper localization capabilities. Zhao et al. [START_REF] Zhao | Robust hashing for image authentication using Zernike moments and local features[END_REF] combined global and local features to construct an image hash. The global features correspond to the Zernike moments of the luminance and chrominance components of the image, while the local ones include the positions and the texture information of salient regions. This algorithm can identify the type of image tamper as well as the modified areas' location. Other robust feature extraction methods for constructing image hashes were also reported including the random Gabor filtering [START_REF] Li | Robust Image Hashing Based on Random Gabor Filtering and Dithered Lattice Vector Quantization[END_REF], the ring partition [START_REF] Tang | Robust image hashing using ring-based entropies[END_REF] and shape contexts [START_REF] Lv | Perceptual Image Hashing Based on Shape Contexts and Local Feature Points[END_REF].

When applied to color images, most of the above schemes convert three color channels (i.e., Red, Green, Blue or RGB) into gray-scale images while discarding the chrominance information. However, exploiting the chrominance information may not only improve the detection performance [START_REF] Zhao | Robust hashing for image authentication using Zernike moments and local features[END_REF][START_REF] Qin | Robust image hashing using non-uniform sampling in discrete Fourier domain[END_REF], but may also make image forgery more difficult due to the fact that the hash contains both luminance and chrominance information. Quaternions can offer a sound way to simultaneously deal with the three color channels without discarding the chrominance information.

They have already been successfully employed in color image registration [START_REF] Wang | Color image registration based on quaternion Fourier transformation[END_REF], image analysis [START_REF] Kantor | Hypercomplex numbers: an elementary introduction to algebras[END_REF][START_REF] Sangwine | Fourier transforms of colour images using quaternion or hypercomplex, numbers[END_REF][START_REF] Ell | Quaternion-Fourier transforms for analysis of two-dimensional linear time-invariant partial differential system[END_REF][START_REF] Pei | Efficient Implementation of Quaternion Fourier Transform, Convolution, and Correlation by 2-D Complex FFT[END_REF][START_REF] Ell | Hypercomplex Fourier transforms of color images[END_REF] and watermarking [START_REF] Tsui | Color image watermarking using multidimensional Fourier transforms[END_REF][START_REF] Lang | A self-adaptive image normalization and quaternion PCA based color image watermarking algorithm[END_REF][START_REF] Chen | Full 4-D quaternion discrete Fourier transform based watermarking for color images[END_REF]. Recently, quaternion discrete Fourier transform (QDFT) was used to generate image hashing and applied to image retrieval [START_REF] Laradji | Perceptual hashing of color images using hypercomplex representations[END_REF]. By following the same path, we proposed a novel image hashing method based on QDFT and log-polar transform for image authentication. QDFT can handle simultaneously the three channels (RGB) of color image without discarding chrominance information. Similar to DFT, the low-frequency coefficients of QDFT contain the main energy of the image and represent essential image features.

In addition, QDFT can be combined with the log-polar transform so as to achieve a set of features that are rotation invariant. We propose thus to take advantage of QDFT used here to build a novel and compact image hash that is robust to content-preserving operations, geometric attacks while being sensitive to malicious tampering operations.

The rest of this paper is organized as follows. Section 2 gives a brief overview of quaternion and QDFT. Section 3 introduces the proposed scheme including its pre-processing, image feature extraction, and hash construction steps as well as its image authentication procedure. Experiments and comparison results are provided in Section 4. Some concluding remarks are given in Section 5.

PRELIMINARIES

Quaternion

Quaternion is a generalization of complex numbers and was introduced by Hamilton in 1843 [START_REF] Kantor | Hypercomplex numbers: an elementary introduction to algebras[END_REF]. A quaternion number has four parts: one real part and three imaginary parts, and can be written as follows:

q a b c d = + + + i j k , (1)
where a , b , c and d are real numbers, i, j and k are imaginary units obeying the following rules:

2 2 2 1 = = = - i j k , , , i = -= = -= = -= j ji k jk kj i ki ik j . ( 2
)
From Eq. ( 2), it can be found that the multiplication rule of quaternion is not commutative. The conjugate and modulus of a quaternion q are respectively defined as follows:

q a b c d = --- i j k , 2 2 2 2 q a b c d = + + + . (3)
A quaternion q with a zero real part is called a pure quaternion. One pixel of a color image f at the spatial position (x, y) has three components, and can be represented as a pure quaternion [START_REF] Sangwine | Fourier transforms of colour images using quaternion or hypercomplex, numbers[END_REF]:

( , ) ( , ) ( , ) ( , ) R G B f x y f x y f x y f x y = + + i j k , (4)
where ( , ) R f x y , ( , ) G f x y and ( , )

B f x y are respectively the RGB component values of f (x, y).

Quaternion discrete Fourier transform

Quaternion or hypercomplex Fourier transform was first introduced in the image processing community by Ell [START_REF] Ell | Quaternion-Fourier transforms for analysis of two-dimensional linear time-invariant partial differential system[END_REF]. Due to the non-commutative property of the quaternion multiplication, there are three different types of QDFT, namely, right-side, left-side and two-side [START_REF] Pei | Efficient Implementation of Quaternion Fourier Transform, Convolution, and Correlation by 2-D Complex FFT[END_REF]. Since that the right side QDFT can be processed in a similar way as the left side QDFT, and the operations are much easier with the left side DQFT or right side QDFT than that with the two-side QDFT [START_REF] Wang | Color image registration based on quaternion Fourier transformation[END_REF][START_REF] Chen | Full 4-D quaternion discrete Fourier transform based watermarking for color images[END_REF], we decided to use the left-side QDFT. For a color image f(x, y) of size M × M, the QDFT and its inverse transform IQDFT are defined as follows [START_REF] Ell | Hypercomplex Fourier transforms of color images[END_REF]:

( )

1 1 2 0 0 1 ( , )
,

ux vy M M M M x y F u v e f x y M µπ - -- + = = = , (5)
( )

1 1 2 0 0 1 , ( , ) ux vy M M M M u v f x y e F u v M µπ - - + = = = , (6)
where is any unit pure quaternion that can be defined as a linear combination of i, j, and k such as

: = i + j + k, , , R, | | = 1.

PROPOSED SCHEME

Our robust image hashing scheme is shown in Fig. 1. Like any robust hashing strategy, it includes an image pre-processing stage followed by an image feature extraction process and a hash construction. We detail each of these stages below.

Pre-processing

The pre-processing procedure we use is illustrated in Fig. 2 for Lena test image. This one takes I 0 as an input image which is first rescaled to a fixed size M×M. This step ensures that the generated image hash will be of fixed length while making it robust against image scaling operation. In this paper, the value M = 256 is chosen. In a second step, to make the generated image hash more robust, a smoothing (i.e. averaging filter) operation based on a k×k window is applied. Its aim is to preserve the essential structures while removing insignificant details without endangering the image content (see Fig. 2(b)). Then, considering that an image can be rotated, we propose to only consider the smoothed pixels within the image inscribed circle as illustrated in Fig. 2 (d). This is obtained by setting to zero the pixels outside the inscribed circle. As exposed in Fig. 2 (d

Feature extraction

Robustness to image scaling being handled in the pre-processing stage, the image features we are looking for have to be robust to image rotation. To do so, we propose to first transform the pre-processing image into log-polar domain and then to apply the QDFT.

Let f 1 (x, y) be a rotated version of the original color image f 0 (x, y) with rotation angle 0 :

1 0 0 0 0 0 ( , ) [( cos sin ), ( sin cos )] f x y f x y x y θ θ θ θ = - + . ( 7
)
If we use the log-polar coordinates: the polar angle in the log-polar system, e the base of the log, the origin (x 0 , y 0 ) being located at the center of the image. Then, Eq. ( 7) can be rewritten as [START_REF] Wang | Color image registration based on quaternion Fourier transformation[END_REF]: respectively. As it can be seen, they are similar except for a cyclical translation along the horizontal axis. When both sides of Eq. ( 9) are transformed by means of QDFT using Eq. ( 5), we obtain:

1 0 0 ( , ) [ , ( )] f f ρ θ ρ θ θ = + . ( 9
0 2 1 0 ( , ) ( , ) v M F u v e F u v θ µπ - = , ( 10
)
where F 0 (u, v) and F 1 (u, v) respectively denote the frequency coefficients of QDFT before and after image rotation, (u, v) denotes the frequency coordinates, and |.| represents the magnitude operator obtained by means of (3), 0 is rotation angle, and M denotes the size of image.

Eq. [START_REF] Swaminathan | Robust and secure image hashing[END_REF] indicates that 1 0 ( , ) ( , )

F u v F u v = .
In other words, the QDFT magnitudes of both f 1 ( , ) and its translated versions f 0 ( , + 0 ) are invariant to rotation. Due to the fact that the low-frequency components of QDFT hold the major part of the image energy and represent the main information on the image content, like for the DFT, we propose to use the magnitude of these low-frequency coefficients as features. More precisely, p magnitude coefficients of a square region from the low-frequency spectrum are selected for that purpose. For example, after the low-frequency coefficients are shifted into the center, if the center of QDFT coefficients is located at u = 129 and v = 129, and p is set to the size 5×5, then the range of selected low-frequency coefficients is:

u, v ∈ [127,131].
As stated in Ref [START_REF] Wang | Attacks on perceptual image hashing[END_REF], if the feature extraction part of a hashing scheme is not dependent of the knowledge of a secret key, or if it is not secured, an adversary will be able to construct a visually different image while reserving the original hash, or able to destroy the hash by introducing small distortion into the image, i.e, preserving its content. In order to secure our perceptual hash extraction, we apply the Arnold transform [START_REF] Chen | Full 4-D quaternion discrete Fourier transform based watermarking for color images[END_REF][START_REF] Arnold | Ergodic problem of classical mechanics[END_REF], parameterized by the secret key K, so as to secretly scrambled the previously p selected QDFT coefficients.

To sum-up, both QDFT and log-polar transform are used in our scheme to capture image features onto which our image hash is built.

Hash construction

Let us represent the p QDFT magnitude coefficients of the selected and scrambled square region into a row vector

(1) (2) ( ) [ , ,..., ] p Z z z z =
, where z (i) corresponds to the i th QDFT magnitude coefficient of vector Z, and i =1, …, p. Each z (i) is determined from the square region according to the selection sequence from top to bottom and from left to right. The robust hash H G is then generated as follows [START_REF] Qin | Robust image hashing using non-uniform sampling in discrete Fourier domain[END_REF]:

( ) ( 1) ( ) ( 1)
1, 0, ( ) 1,..., -1. 0, 0,

i i G i i if z z H i i p if z z + + - ≥ = = - < (11)
In order to enhance the security of H G , this one is secretly scrambled based on a process parameterized by a secret key K H . Without the knowledge of K H , one will not be able to generate a valid hash. Here, we used the MatLab function randperm () to encrypt. In fact, it is a randomly permutation problem. According to the secret key, it generates a scrambled serial number. The encrypted hash is formed by the scrambled serial number.

Image authentication

In a classical scenario, the original image I 0 is transmitted along with its robust hash H G0 to a recipient through the Internet or a third party certification organization. At the reception, the recipient can then check the image authenticity. The received image 1 I , a version of I 0 that may have undergone some content-preserving or malicious tampering operations, is utilized to generate the hash 1 G H according to the same procedure described above. These two hashes, 0 G H and 1 G H , are then compared to determine whether the received image is tampered or not. The normalized Hamming distance D 1 is used in our system to measure the similarity between the two hashes:

1 1 1 0 1 0 1 1 1 ( , ) ( ( ) ( ( ))) 1 p G G G G i D d H H Key H i Key H i p - = = = ⊕ - ( ) , (12)
where ' ⊕ ' represents the exclusive-or operation, and Key(.) denotes the encryption operation. If the distance D 1 is less than a pre-fixed threshold , the two images are considered as similar images; otherwise, they are viewed as distinct images or tampered images. The determination of the parameters of our robust image hashing scheme is analyzed in the experimental section.

EXPERIMENTAL RESULTS

Image data sets, performance criteria

To evaluate the performance of the proposed scheme, the UCID [START_REF] Schaefer | UCID -An Uncompressed Colour Image Database[END_REF] image database (including 1338 color images of size 512×384 or 384×512 pixels) was selected. In all the experiments, the first 1000 images of this database were used unless stated otherwise. Some example images are displayed in Fig. 4.

The receiver operating characteristics (ROC) [START_REF] Fawcett | An introduction to ROC analysis[END_REF] is used to evaluate the influence of the parameters and the authentication performance in comparison with other methods. The true positive rate (P TPR ) and the false positive rate (P FPR ) of the ROC curve indicate the perceptual robustness and the discriminative capability (or anti-collision capability) of an image hash, respectively. They are defined as follows:

1 1 1 ( ) ( ) TPR n D P N λ λ < = , 2 1 2 ( ) ( )
FPR n D P N λ λ < = (13)
where n 1 corresponds to the amount of pairs of visually identical images classified into similar images, n 2 is the amount of pairs of different images classified into similar images, while N 1 and N 2 denote the total number of pairs of visually identical and different images, respectively. In this work, we experimentally identify the values of these parameters so as to ensure a good trade-off in terms of detection performance of our system. To do so, the first 1000 of the UCID images were used in this experiment as original image database, each image being submitted to 58 content-preserving operations, the parameterization of which is given in Table 1. These operations include pepper & salt noise, Gaussian noise, brightness adjustment, different filters, JPEG compression, cropping, scaling and rotation. There are thus a total of N 1 =58000 visually identical images pairs, which are looked as constituting the visually identical image database. The normalized Hamming distance D 1 between the original image database and the visually identical image database is computed [START_REF] Zhao | Robust hashing for image authentication using Zernike moments and local features[END_REF]. Note that there are N 2 = 1000× (1000-1)/2 =499500 pairs of different images.

(1) Influence of the pre-processing average filter window size (k). To evaluate the influence of the pre-processing stage average filter and to select its suitable window size, we calculated the average Hamming distance between hashes on the whole image test set (i.e. 499500 pairs of different images and of 58000 pairs of similar images) for different window sizes k = 3, 5, 7, and 11 while considering a larger number of QDFT coefficients p (fixed to p=1000). These distances are shown in Table 2, where DD1 and DD2 denote the mean of normalized Hamming distance between hashes of visually similar images and distinct images, respectively. Small values of DD1 reflect a good robustness while large values of DD2 mean a good discriminative capability. We also introduce DD3 which corresponds to the absolute difference between DD2 and DD1. A large value of DD3 indicates that the corresponding average filter window size achieves a good trade-off between robustness and discriminative capability.

As it can be seen, if the averaging filter window is of very small size, e.g. 3×3, the value of DD1 is large (see bold numbers in Table 2) making our hash non-robust to content preserving operations like average and median filtering. When increasing the size of the filter window, DD1 values decrease whatever the image modification. The averaging filter induces an image information loss which influences the robustness of our hash. At the same time, this filtering operation also impacts the discrimination capability of our scheme as shown by DD2 values. These values decrease when increasing the window size. However, a good tradeoff in terms of robustness and discrimination capability of our hash can be achieved by looking at DD3, the highest value of which DD3=0.373 is obtained for a window of 7×7 pixels. Based on these experiments, we decided to fix the averaging filter window size to k = 7.

(2) Influence of the number of selected low-frequency QDFT coefficients (p). As exposed in Section 3.3, the length of our robust hash depends on the parameter p. In general, the smaller the value p is, the better the perceptual robustness will be but the worse the anti-collision performance is, and vice-versa. In order to find an appropriate trade-off, the previous test color images we built from the UCID ones were authenticated considering different values of p within the set 5 2 , 9 2 , 13 2 , 17 2 , 21 2 , 25 2 and 51 2 . Fig. 5 shows the corresponding ROC curves, with a zoom in Fig. 5(b) of Fig. 5(a). Slight improvements are observed when p is increasing from 81 to 289. For higher values the performance decreases. Thus, one can see that a good trade-off between the perceptual robustness and the anti-collision capability of our hash is obtained for p values in the range [169,289]. In the rest of this paper, the value p=15 2 =225 was retained in order to have an appropriate short hash length. Indeed, in that case, H G is of 224 bits according to Eq. ( 11). distances for different image pairs are larger than 0.2, while almost all distances between identical images are less than 0.2. This is the reason why the authentication threshold we use in the sequel is set to = 0.2.

Performance test

As described in the introduction, robust image hashing should contain some key properties: perceptual robustness, anti-collision capability, sensitivity and security. In this section, we examined these aspects so as to establish the performance of our hash.

(1) Perceptual robustness. Perceptual robustness means that both the received image and the original image have a very high probability of generating similar hash. Fig. 6(a) shows the perceptual robustness of the proposed method considering the above parameterization. It can be seen that only a few image pairs (about 0.2%) are looked as different images when is set to 0.2. Notice that these misclassified images are due to strong parameter attacks such as pepper and salt noise with density 0.1 and a median filter with a window size 11×11. We further use the well-known images such as Lena, Airplane, Pepper, Mandrill and Sailboat and lake to show the robustness.

Almost all attacks described in Table 1 were applied. The normalized Hamming distance values corresponding to each operation are displayed in Fig. 7. It can be seen that the normalized Hamming distance DD1 is less than 0.2.

This perceptual robustness can be achieved for three main reasons: (i) the generated hash utilizes the low-frequency QDFT coefficients whose interrelation is hardly changed by content-preserving operations; (ii) the average preprocessing filter makes the method more robust against content-preserving operations such as noise, filtering operation and JPEG compression; (iii) the magnitude of QDFT coefficients of the log-polar transformed image is invariant to rotation. The proposed method is thus perceptually robust not only to geometric attacks but also to content-preserving operations. The index of wrong secret key Normalized hamming distance Fig. 8. Distribution of normalized Hamming distance when using 50000 wrong secret keys (3) Security. To validate the security of the system, we tested 50000 wrong secret keys. The hash distance distribution of image pairs between right and wrong secret keys K H is shown in Fig. 8. We observe that all these distances are higher than 0.2. If one does not know the secret key and the hash algorithm, the probability of generating the same hash is about 1/2 224 since our hash length is 224 bits. If one does not know the secret key but know the hash algorithm, it is a special permutation problem of 224 bits. Therefore, the proposed method meets also the security requirements.

(4) Sensitivity. When some tampering operations are applied to an image, the corresponding image hash should be distinct from the original image one. To validate the sensitivity of the proposed method to such visual content alterations, we have simulated two kinds of tampering operations. They rely on copying a portion of a "source image" and pasting it into a "destination image" which can be same as the source image ("copy-move attack") or another image ("cut-and-paste attack"). As illustrated in the second and third row of Fig. 9, the upper left corner of the pasted area corresponds to the destination image center. We have considered different sizes of the copied area so as it correspond to 4%, 9%, 16% and 25% of the destination image dimensions. It can be seen from Fig. 9 that for such kind of attack the resulting normalized Hamming distances DD2 are higher than 0.2. As shown, our method is sensitive to content alterations taking advantage of the fact that the magnitude of QDFT coefficients contains both the luminance and the chrominance information. To further validate the sensitivity of our hash to visual content alterations, the first 1000 images from the UCID database were also tested. Similarly to the above copy-move and cut-and-paste attacks, tampering ratios that were considered are of 4% 9% 16% 25% 36% and 49%, leading thus to 6000 normalized Hamming distances shown in Fig. 10. We observe that all distances are larger than 0.2 except for one image where the tampering ratio is equal to 4%. This confirms that the proposed method is sensitive to content alterations when the tampering ratio is higher than 4%. In addition, we have also conducted experiments with a tampering position near the boundary (a tampering ratio about 10% was chosen) as shown in the fourth row of Fig. 9. We observe that such operation leads to a sensitivity decrease (DD2 is higher than 0.2 but close). The number of forged images Normalized hamming distance Forged ratio 4% Forged ratio 9% Forged ratio 16% Forged ratio 25% Forged ratio 36% Forged ratio 49% In order to better evaluate the performance of our QDFT based scheme, we compared it with different methods based on DFT [START_REF] Qin | Robust image hashing using non-uniform sampling in discrete Fourier domain[END_REF], QDFT [START_REF] Laradji | Perceptual hashing of color images using hypercomplex representations[END_REF], Zernike moments [START_REF] Zhao | Robust hashing for image authentication using Zernike moments and local features[END_REF] and NMF [START_REF] Tang | Robust perceptual image hashing based on ring partition and NMF[END_REF]. On one hand, these methods represent the current state of the art and, on the other hand, their authors suggested working on color images. In order to further examine the influence of the log-polar transform on the performance of scheme; three QDFT based schemes were compared. The first one, denoted by QDFT1, corresponds to the method already described. The second one, denoted by QDFT2, is identical to QDFT1 except that the log-polar transform is omitted. The parameters used for these two methods are the same. The QDFT method reported in [START_REF] Laradji | Perceptual hashing of color images using hypercomplex representations[END_REF] is referred to as QDFT3. As previously, the first 1000 images of the UCID image database were used to generate similar image pairs and distinct image pairs.

The comparison is worked out using ROC curve, which is a good measure for a fair comparison since it allows comparing methods based on different distance measurement criterion (i.e. Zernike moments [START_REF] Zhao | Robust hashing for image authentication using Zernike moments and local features[END_REF] use Euclidean distance, NMF [START_REF] Tang | Robust perceptual image hashing based on ring partition and NMF[END_REF], the correlation coefficients, DFT [START_REF] Qin | Robust image hashing using non-uniform sampling in discrete Fourier domain[END_REF], QDFT [START_REF] Laradji | Perceptual hashing of color images using hypercomplex representations[END_REF] and our scheme, the Normalized Hamming distance). The ROC curves are displayed in Fig. 12 (a). It can be seen that both QDFT1 and QDFT2 leads to a better behavior than the other hashing methods, QDFT3 shows the worst performance. The main reason stands on the way this hash is constructed. Basically, it directly divides the frequency domain into block of size 8×8, and then generates the hash component value by comparing the average energy of block with the one of the whole image.

But due to the fact about 9% of low-frequency coefficients carry more than 99.5% of the image energy as stated in [START_REF] Wang | Attacks on perceptual image hashing[END_REF], most of hash component values are '0', and only a few ones have value '1'. This leads to a very bad anti-collision capability and also impacts the overall performance of this scheme.

Fig. 1 .

 1 Fig. 1. Block diagram of the proposed robust image hashing scheme.

  ) and Fig.2(e), the image and its rotated version have the same pixel values. This allows us minimizing the influence of information loss, especially the loss of pixels in the image corners, in the case of image rotations with large angles.

Fig. 2 .

 2 Fig. 2. Pre-processing procedure applied to an image and to a rotated version of an angle of 45°. (a) Image I 0 , (b)-(c) Smoothed images, (d)-(e) inscribed circle in smoothed images (b) and (c).

  the logarithm of the radial distance from the origin (x 0 , y 0 ) of the image,

Fig. 3 .

 3 Fig. 3. Results of Log-Polar Transform applied to Fig. 2(d) and (e)

Fig. 4 .

 4 Fig. 4. Examples of test image

  ROC curve of different p values (b) Zoom-in view of the ROC curves Fig. 5. ROC curves for different values p .

6 .( 3 )

 63 Distribution of the normalized Hamming distance 1 D Authentication threshold . Fig. 6(a) and (b) show respectively the Hamming distance distributions for hashes of visually identical color images and for hashes of visually different images. It can be seen that almost all

Fig. 7 .

 7 Fig. 7. The distribution of distances DD1 in different content-preserving operations

Fig. 9 .

 9 Fig. 9. Examples of original and tampered images. Top: original images. Second row: tampered images after a copy-move attack of size 4%, 9%, 16%, 25%, obtained normalized Hamming distance are DD2=0.3259, 0.3527, 0.3839, 0.4107, respectively. Third row: tampered image after a cut-and-paste attack of size 4%, 9%, 16%, 25%, obtained normalized Hamming distance values are DD2=0.2902, 0.4107, 0.4107, 0.4241, respectively. Fourth row: tampered image after a tampering rate of 10%. The resulting normalized Hamming distance values are DD2=0.2679, 0.2143, 0.2455, 0.3438, respectively.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Distribution of normalized Hamming distances for 6000 different tampering image pairs

  Gaussian filtering (g) Median filtering

Acknowledgments

This work was supported by the National Basic Research Program of China under Grant 2011CB707904, the National Natural Science Foundation of China under Grants 61073138, 61103141, 61271312, 61201344, the Research Fund for the Doctoral Program of Higher Education of Ministry of Education of China under Grants 20110092110023 and 2012009212003, the Qing Lan Project, DZXX-031, BY2014127-11 of Jiangsu province, and the Natural Science Foundation of the Jiangsu Higher Education Institutions of China under Grant 13KJB520015.

From these results, it appears that these methods are robust to some content-preserving attacks such as Gaussian noise, Gaussian filtering, average filtering, JPEG compressing, etc. However, our method behaves better on the overall. For instance, Fig. 12(d)-(e) shows that the Zernike moment method is sensitive to brightness and contrast adjustment. Fig. 12(e), (g), (j) and (k) show that the QDFT3 is not robust to contrast adjustment, median filtering, rotation and cropping attack. For the rotation operations illustrated in Fig. 12(j), the NMF method and QDFT1 method obtain better performance than the other four methods. When comparing QDFT1 and QDFT2, the former shows better robustness to rotation and cropping attacks since it employs the log-polar transform in feature extraction process. Thus, the proposed method is robust to almost all content-preserving operations. This is due to the joint combination of (i) the low-frequency magnitude coefficients of QDFT, (ii) the log-polar transform of the image central part, and (iii) the filtering operation enhances its robustness and increases its discriminative capability.

Sensitivity comparison

To do that, we have built three test image sets including three tampering situations: copy-move, cut-and-pasting and region color changing. Since it is a tedious and time-consuming work to generate the tampering images using Photoshop soft, the copy-move and cut-and-pasting tampering images were automatically generated by computer.

The first test set (Data 1) contains 2000 tampering images as used in subsection 4.3 for copy-move and cut-and-pasting tampering operations. The second test set (Data 2) differs from the first test set by using the 10% ratio of tampering areas with a randomly selected position of pasting. The third test set (Data 3) contains 150 tampering images. 50 region color changing tampering images were generated by Photoshop, another 100 tampering images were selected from other two data sets (each one including 50 images).

Fig. 13 summarizes the results obtained. For Data 1, the proposed method provides the best rate of recognition with a better sensitivity to slight tampering area (4%) than the other three methods. For Data 2, the DFT method [START_REF] Tang | Robust perceptual image hashing based on ring partition and NMF[END_REF] leads to the best results. It remains insensitive to tampering position changes when the tampering area is larger than 10%. For Data 3, the Zhao's method [START_REF] Zhao | Robust hashing for image authentication using Zernike moments and local features[END_REF] provides the best rate of recognition in the case of region color changing.

Data1

Time complexity and qualitative comparison

To synthetize the comparison of our method and other hashing methods, a qualitative overview with the current state of art is given in Table 3. We observe that the hash of our method has the shortest length and that it is robust against various content-preserving operations and any angle rotation operations. We also consider the average time required for generating a hash value in our experimental environment (i.e. CPU 3.2 GHz, 4G memory and running MatLab 7.12). Since the local features of Zhao's method [START_REF] Zhao | Robust hashing for image authentication using Zernike moments and local features[END_REF] are mainly used to classify and locate a tampered region, the time computation of the global features generating the hash is only considered. It can be seen that the time running for the five methods is very short and less than 0.22 s. The DFT method [START_REF] Tang | Robust perceptual image hashing based on ring partition and NMF[END_REF] is the fastest and needs only 0.07 s.

CONCLUSION

In this paper, we have developed a novel robust image hashing method devoted to color images using QDFT and log-polar transforms. By exploiting the QDFT into the log-polar domain, we better take into account the three image color planes and the resulting hash is robust to rotation attacks as well as to common image content-preserving operations. The parameters' values of our approach have been determined by means of intensive experiments on a large image data set so as to establish a good tradeoff in terms of hash robustness against content-preserving operations (JPEG compression, brightness adjustment, average and median filtering and large angle rotation attacks) and a good sensitivity to non-authorized image content alterations (copy-move and copy-paste attacks). Compared to the current state of art, for a hash length of 224 bits our scheme provides on the overall a better performance. Our future work will focus on designing a robust hash capable to locate tampered regions of an image as well as to determine the type of the tampering.


Mahmoud Hassan
email: mahmoud.hassan@univ-rennes1.fr

Pascal Benquet

Arnaud Biraben

Claude Berrou

Olivier Dufor

Fabrice Wendling

Dynamic reorganization of functional brain networks during picture naming

For efficient information processing during cognitive activity, functional brain networks have to rapidly and dynamically reorganize on a sub-second time scale. Tracking the spatiotemporal dynamics of large scale networks over this short time duration is a very challenging issue. Here, we tackle this problem by using dense electroencephalography (EEG) recorded during a picture naming task. We found that (i) the picture naming task can be divided into six brain network states (BNSs) characterized by significantly high synchronization of gamma (30-45 Hz) oscillations, (ii) fast transitions occur between these BNSs that last from 30 ms to 160 ms, (iii) based on the state of the art of the picture naming task, we consider that the spatial location of their nodes and edges, as well as the timing of transitions, indicate that each network can be associated with one or several specific function (from visual processing to articulation) and (iv) the comparison with previously-used approach aimed at localizing the sources showed that the network-based approach reveals networks that are more specific to the performed task. We speculate that the persistence of several brain regions in successive BNSs participates to fast and efficient information processing in the brain.

Introduction

Any cognitive process involves the activation of a large-scale functional brain network [START_REF] Bressler | Large-scale brain networks in cognition: emerging methods and principles[END_REF]. In visual, attentional and memory processes, this network is characterized by increased synchronization of cortical oscillations (in the gamma frequency range [START_REF] Doesburg | Large-scale gamma-band phase synchronization and selective attention[END_REF], in particular but not only) across distant neuronal assemblies distributed over distinct brain areas.

The accurate tracking of the spatiotemporal dynamics of large-scale networks over the duration (often as short as a few hundreds of ms) of cognitive processes is still a challenging issue [START_REF] Allen | Tracking wholebrain connectivity dynamics in the resting state[END_REF][START_REF] Hutchison | Dynamic functional connectivity: promise, issues, and interpretations[END_REF]. A number of theories have been elaborated to explain these spatiotemporal dynamics. It has been hypothesized that functional brain networks engage in fast transitions between transiently stable states, each characterized by a network with intrinsic dynamics and with specific functional relationships between neuronal assemblies [START_REF] Hansen | Functional connectivity dynamics: Modeling the switching behavior of the resting state[END_REF][START_REF] Sporns | Networks of the Brain[END_REF]. According to this theory, the substrate of cognitive processes would correspond to a sequence of switches between networks and, thus, to time-and space-dependent fluctuations in the node and edge properties of the global network.

The validation of such hypotheses for task-related data requires the following of brain processes at the millisecond time-scale. This can barely be achieved using fMRI data for a simple and well-known reason: although they are characterized by an excellent spatial resolution, BOLD signals reflect the metabolic and hemodynamic response of neuronal assemblies (at voxel level). This slow response (seconds) is obviously related to the fast dynamics of cortical oscillations taking place over interconnected neuronal assemblies and defining functional networks, but indirectly (i.e. through the neuro-glial-vascular coupling [START_REF] Logothetis | Neurophysiological investigation of the basis of the fMRI signal[END_REF]).

In this study, we address this issue using EEG source connectivity analysis to track the spatiotemporal dynamics of large-scale networks associated with cognitive activity. We collected dense-EEG data from 21 subjects performing a picture naming task. We then reconstructed the functional networks in both their spatial and temporal dimension, over the entire duration of the cognitive process (from image perception to motor response) using a recently reported method [START_REF] Hassan | EEG Source Connectivity Analysis: From Dense Array Recordings to Brain Networks[END_REF]) that combines i) the solution to the inverse EEG problem, ii) the estimation of brain connectivity from phase locking values and iii) the segmentation of functional networks using a clustering method [START_REF] Mheich | A new algorithm for spatiotemporal analysis of brain functional connectivity[END_REF] (see figure 1A). Our results reveal that the cognitive process can be decomposed into a sequence of transiently-stable and partiallyoverlapping networks. We assume, based on the state of the art of the picture naming task, that each network might be associated with a specific function [START_REF] Levelt | An MEG study of picture naming[END_REF]) (visual percept computing, lexical concept activation, selecting the target word from the mental lexicon, phonological encoding, phonetic encoding, and initiation of articulation) of the whole cognitive process. The results show that dense-EEG can bring highly valuable information about cortical networks, with both high spatial (1000 cortical regions) and temporal (ms time-scale) resolution. We speculate that the identified brain network states (BNSs) contribute to fast and efficient information processing in the brain.

Materials and methods

Picture naming task

Twenty one right-handed healthy volunteers (11 women: mean age 28 year; min: 19, max: 40 and 10 men: mean age 23 years; min: 19, max: 33), with no neurological disease, were involved in this study. Participants were asked to name at a normal speed 148 displayed pictures on a screen using E-Prime 2.0 software (Psychology Software Tools, Pittsburgh, PA) [START_REF] Schneider | E-Prime: User's guide: Psychology Software Incorporated[END_REF]. The images were selected from a database of 400 pictures standardized for French [START_REF] Alario | A set of 400 pictures standardized for French: Norms for name agreement, image agreement, familiarity, visual complexity, image variability, and age of acquisition[END_REF] and were used during session about eight minute. They were controlled according to several parameters (see table S1). All pictures were shown as black drawings on a white background. Order of presentation was randomized across participants. Naming latencies were determined as the time between picture onset and the beginning of vocalization recorded by the system. Oral responses were recorded and then analyzed with Praat software to set the voice onset time [START_REF] Boersma | Praat, a system for doing phonetics by computer[END_REF].

This study was approved by the National Ethics Committee for the Protection of Persons (CPP), conneXion study, agreement number (2012-A01227-36), and promoter: Rennes University Hospital. All participants provide their written informed consent to participate in this study. The ethics committee has approved the consent procedure. A typical trial started with the appearance of an image during 3 seconds followed by a jittered inter-stimulus interval of 2 or 3 seconds randomly. Most responses were given while the image was still present on the screen. Errors in naming were discarded for the subsequent analysis. A total of 2926 on 3108 events were considered. The fastest response time delay for an event was < 600 ms (see figure S1 for the responses delays of all subjects).

Data

The brain activity was recorded using dense-EEG, 256 electrodes, system (EGI, Electrical Geodesic Inc.). The main feature of this system is the large coverage of the subject's head by surface electrodes allowing for the improved analysis of the intracerebral activity from noninvasive scalp measurements, as compared with 32 -to 128-electrodes standard systems. EEG signals were collected with a 1 kHz sampling frequency and band-pass filtered between 3 and 45Hz. Each trial was visually inspected, and epochs contaminated by eye blinking, movements or any other noise source were rejected and excluded from the analysis performed using the EEGLAB open source toolbox [START_REF] Delorme | EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis[END_REF].

EEG source connectivity

A crucial step when realizing EEG source connectivity analysis is the choice of three factors: the method used to solve the inverse problem, the method used to compute the functional connectivity among the time series of the reconstructed sources and the number of electrodes used on the scalp. Very recently, we have described a comparative study of these factors and we showed that a combination of the weighted Minimum Norm Estimate (wMNE) with the Phase Locking Value (PLV) using high resolution EEG is the best combination among the tested combination [START_REF] Hassan | EEG Source Connectivity Analysis: From Dense Array Recordings to Brain Networks[END_REF]. This combination was used in the presented work. where G and B(t) are respectively the matrix containing the lead fields of the dipolar sources and the additive noise. In the general case, the inverse problem consists in finding an estimate D (t) of the dipolar source parameters (typically, the position, orientation and magnitude), given the EEG signals S(t) and given the gain matrix G. This matrix can be computed from a multiple layer head model (volume conductor) and from the position of electrodes. For instance, the Boundary Element Method is a numerical method classically used in the case of realistic head models.

As this problem is ill-posed (P>>Q), physical and mathematical constraints have to be added to obtain a unique solution among the many solutions that minimize the residual term in the fitting of measured EEG signals. Using segmented MRI data, the source distribution can be constrained to a field of current dipoles homogeneously distributed over the cortex [START_REF] Dale | Improved Localization of Cortical Activity by Combining EEG and MEG with MRI Cortical Surface Reconstruction: A Linear Approach[END_REF], and normal to the cortical surface.

Technically, in the source model, we assumed that EEG signals are generated by macrocolumns of pyramidal cells lying in the cortical mantle and aligned orthogonally with respect to its surface [START_REF] Nunez | Electric fields of the brain: the neurophysics of EEG[END_REF]. Thus, the electrical contribution of each macro-column to scalp electrodes can be represented by a current dipole located at the center of gravity of each triangle of the 3D mesh and oriented normally to the triangle surface. Using this source space, the weighted Minimum Norm Estimate (wMNE) method only estimates the moment of dipole sources. The wMNE compensates for the tendency of classical MNE to favor weak and surface sources. This is done by introducing a weighting matrix W S :

T 1 T wMNE S S D G W G I G W S    ˆ()
where matrix W S adjusts the properties of the solution by reducing the bias inherent to MNE solutions. Classically, W S is a diagonal matrix built from matrix G with non-zero terms inversely proportional to the norm of the lead field vectors. The value of  is computed relatively to the signal to noise ratio for each signal computed as the ration between the poststimuli period to the pre-stimulus (200 ms). The  value was between 0.1 and 0.3.

The sources were reconstructed for each trial (same number of sources for each trial) and the functional connectivity is then computed between the reconstructed sources using the phase synchronization (PS) method. The first step for estimating the PS is to extract the instantaneous phase of each signal. We used the method based on Hilbert transform in our study. The second step is the definition of an appropriate index to measure the degree of synchronization between estimated instantaneous phases. To measure PS, the phase locking value (PLV) method is used as described in (Lachaux, Rodriguez, Martinerie, & Varela, 1999). For each source pair, x and y, at time t (t = t 1 ,. ., t T where T= D * f s ; D and f s denote the signal length relative to the onset and the sampling frequency, respectively) for the Tr trials and for subject j (j = 1...M, where M denotes the number of subjects), PLV is defined as: . The functional connectivity was computed at the low gamma band . This frequency band is the most relevant one in the context of the similar cognitive task performed by the subjects, as reported in [START_REF] Fell | The role of phase synchronization in memory processes[END_REF][START_REF] Rodriguez | Perception's shadow: long-distance synchronization of human brain activity[END_REF][START_REF] Supp | Directed cortical information flow during human object recognition: analyzing induced EEG gamma-band responses in brain's source space[END_REF]. The PLVs were then averaged over subjects:

N N t t T  , computed
for the V pairs of x and y sources, where V is equal to (Nc.(Nc-1)/2)) and Nc is the number of ROIs.

Functional connectivity states

We recently developed an algorithm to decompose cognitive task into functional connectivity states [START_REF] Mheich | A new algorithm for spatiotemporal analysis of brain functional connectivity[END_REF]. The objective of this algorithm is to identify clusters among the T networks () Nt . The proposed algorithm is based on the K-means clustering of the connectivity networks obtained by the PLV method. This approach allowed us to summarize brain networks into a limited number of dominant networks over given time period.

Briefly, the averaged connectivity matrices over all subjects were first obtained. K networks (varies from 3 to 12) were then randomly selected and the spatial correlations between the K networks and all the T networks were computed. This gives a spatial correlation value for each K at each time instant and for any of the T networks only one of the K networks produces the highest spatial correlation. Finally, the cross-validation criterion, first introduced by (Pascual-Marqui, [START_REF] Pascual-Marqui | Segmentation of brain electrical activity into microstates: model estimation and validation[END_REF], was used to determine the optimal number of networks that explained the best the ongoing cognitive task. To investigate the inter-subject variability, we added an index called 'network presence' which calculate the ratio (in %) of the identified networks/clusters among all the subjects (see Supplementary document for detailed description of the algorithm).

Statistical tests

The segmentation algorithm based on k-means, produces a number of functional connectivity states (clusters). We assume that these clusters reflect changes in the cognitive state. To verify this assumption, we compare our results with an appropriate null model. The identified clusters are compared with those obtained by shuffling the original data using surrogates (procedure detailed in [START_REF] Hlinka | On the danger of detecting network states in white noise[END_REF]). Briefly, we use multivariate Fourier transform surrogates (n sur =100) generated from the original EEG data for all trials. Such surrogates correspond to realizations of linear stationary process with conserved auto-and cross-correlation characteristics. The same processing steps were performed on these surrogates as on the original data: filtering in the gamma band: 30-45Hz, computing the functional connectivity using phase locking values at this frequency band, thresholding the connectivity matrices and segmenting into clusters using k-means algorithm. We then compared the spatial distributions (Sd) and the temporal profiles (Tp) of the original clusters with those obtained using surrogates.

Concerning the Tp, when a cluster is significant (i.e. related to the cognitive process and not to chance), the value of Tp for the surrogate data (Tp surr ) will be different than the one of the original cluster (Tp org ). The underlying null hypothesis is that the obtained clusters are significantly different from the spurious networks states that can be always observed in random data or white noise [START_REF] Hlinka | On the danger of detecting network states in white noise[END_REF]. The null hypothesis is tested by comparing the Tp surr and Tp org using a statistical test. The "Rank test" is used to reject or accept the null hypothesis. Basically, [Tp org ; Tp surr ] is sorted in increasing order and the rank index for Tp org is returned. With a number of surrogates (n_surr =100 for example), if this rank is > 95 and < 5 (significance level at 95%), this means that it lies in the tail of the distribution, and that the null hypothesis can be rejected (two-tailed test) with a significance of p=2*(1/ (n_surr+1)) =0.019. A similar test is used for Sd. For multiple tests and to deal with the family wise error rate (FWER), we used the well-known Bonferroni correction test, which is considered as the simplest and most conservative method to control the FWER problem

Regions of interest and network measures

We used Freesurfer (Fischl, 2012) to register a labeled mesh from an average brain, where each label corresponded to one of 148 anatomical cortical regions [START_REF] Destrieux | Automatic parcellation of human cortical gyri and sulci using standard anatomical nomenclature[END_REF]. This output provided a standardized partition of the cortex into 148 regional areas. Each of these areas was then subdivided into a set of small sub-regions using Brainstorm [START_REF] Tadel | Brainstorm: a user-friendly application for MEG/EEG analysis[END_REF], resulting in 1000 ROIs covering the whole cortex. This segmentation provided us with high resolution connection matrices (see figure 1A). These ROIs were then considered as In order to obtain a sparse, weighted, undirected graph, the adjacency matrices were thresholded. For each matrix, 10000 edges were retained. All corresponding weight values were positive. Then, the strength measure was used to characterize the nodes in the obtained weighted networks. This measure is defined as the sum of all edge weights for each node. For any node i, the strength w i k is defined as:

w i ij jN kW   
where N is the set of all nodes in the graph and ij W is the connection weight between two nodes i and j. This measure was computed using the Brain Connectivity Toolbox (BCT) [START_REF] Rubinov | Complex network measures of brain connectivity: uses and interpretations[END_REF].

Software

We

Results

Functional connectivity states identified from dense-EEG recordings. We have recently developed a new clustering algorithm to detect the brain network states (BNSs) at millisecond scale from dense EEG recordings [START_REF] Mheich | A new algorithm for spatiotemporal analysis of brain functional connectivity[END_REF] (figure 1B). A BNS is defined by a transiently-stable brain network in which nodes are associated with distributed neuronal assemblies and in which edges denote significantly high phase synchronization of their oscillatory activity in the gamma band (30-45 Hz) (see Materials and Methods for details about the algorithm). In order to assess brain dynamics during information processing (from perception to behavioral response), this algorithm was applied to high-resolution group averaged connectivity matrices obtained from dense EEG signals recorded in 21 subjects performing a picture naming task.

Results are shown in figure 2 and they revealed that the cognitive process could be divided into 6 BNSs. To globally characterize the networks for each of these 6 BNSs, we show the number of nodes for five macroscopic regions (occipital (O), parietal (P), temporal (T), frontal (F) and central (C)) expressed in percent (with respect to the total number of nodes).

The first BNS corresponds to the time period ranging from 0 (stimulus onset) to 119 ms. This network mainly involves the occipital pole (L: 62%, R: 64%), the temporal lobe (L: 30% and R: 14%) and the parietal lobe (L: 1% and R: 22%). A second stable state (BNS2) was then observed between 120 and 150 ms where the network is mainly located in the occipital (L: 88%, R: 79%) and temporal (L: 12% and R: 14%) lobes. The third state (BNS3) was identified over the time period 150-190 ms. It was found to be very similar to the previous network in term of spatial location of nodes. However, a main difference relates to appearance, in BNS3, of nodes located in the right frontal lobe (L: 0% and R: 5%). A dominant network (BNS4) was then later identified over the period 190-320 ms. This network is characterized by long-range connections between the occipital (L: 44% and R: 61%) and the frontal lobe (L: 28% and R: 22%). This network was then followed by BNS5 over the time period 320-480 ms in which connections with the pre-central region appear (L: 18% and R: 13%). Finally, the last network BNS6 was obtained over the period 481-535 ms. This functional network mainly involved the temporal lobe (L: 39% and R: 31%) and central region in each hemisphere (L: 26% and R: 36%). In addition, our results showed also strong inter-hemisphere connections mainly between the occipital lobes at BNS1, BNS2, BNS3 and BNS4 and between frontal lobes at BNS4 and BNS5. (See figure S2 for multiple view -left, right, top, front and back-of the obtained networks).

The identified BNSs were compared with clusters obtained from null model (surrogate data).

Typical examples of the results corresponding to the BNS2 and BNS3 are shown in figure 3.

To evaluate the temporal profiles of the identified BNSs, we compared the duration of both networks with those obtained by the 100 surrogates. Figure 3 shows significant differences between the duration of clusters identified from surrogates data and the original ones (p<0.01). To evaluate the spatial distributions of the BNSs, we chose to compare the number of nodes identified in the occipital regions (involved mainly in the visual information at this period) for the original and surrogates data. Results also show significant differences (p<0.01)

for both BNS2 and BNS3. The same results were obtained for all the other BNSs.

EEG source localization vs. functional connectivity. We compared the above-described source connectivity results with those obtained from the more classical source localization analysis, as reported in Figure 4. The distributions of sources were obtained from dense-EEG recordings using the weighted Minimum Norm Estimate (wMNE) algorithm. They were averaged over subjects and over time, for each of the 6 periods identified with the segmentation algorithm. Results show that both approaches lead to dramatically different results regarding involved brain regions. The main differences are described hereafter. During time periods t1, t2 and t3, the source localization method disclosed a pronounced bilateral activation of the dorsolateral frontal cortex. Interestingly, this region was not present in the networks associated with the BNSs identified over the same three periods by the source connectivity method. To some extent, the same result holds for t4, for the right hemisphere.

Another discrepancy was observed for the inferior occipital and the occipito-temporal regions that were not disclosed by the wMNE algorithm over period t5 while the same regions were not part of the network associated with BNS5 revealed by the source connectivity method.

Finally, over t6, differences between both methods relate to the central and insula regions, as they did not appear in localized sources but were present in the network associated with BNS6.

Spatiotemporal dynamics of brain networks and associated neocortical areas.

In order to further precise characterization of the spatiotemporal dynamics of identified networks and in order to accurately identify the neocortical areas involved in these networks, we combined a quantitative network measure with the high-resolution anatomical parcellation of the brain surface (~ 1000 ROIs). First, the "strength" was used as network measure (high "strength" means "strong" connections at a given node). Second, the parcellation/labelling step allowed us to establish a direct link between spatially distributed nodes showing significantly high strength values (> mean + 1SD), on the one hand, and the corresponding cortical areas (labelled ROIs) on the other hand. Results are presented in figure 5A for the 6 periods identified during the cognitive process. First, it was found that 3 (minimum) to 5 (maximum) brain areas are involved in a given BNS. Second, the highest strength values were observed for t2 (L: 0.85±0.1 and R: 0.71±0.1) and t3 (L: 0.9±0.13 and R: 0.91 ± 0.12) as compared with t1 (L: 0.3±0.08 and R: 0.32 ±0.11), t4 (L: 0.42±0.04 and R: 0.4±0.4), t5 (L: 0.21±0.03 and R: 0.17±0.04) and t6 (L: 0.39±0.08 and R: 0.41±0.071). Third, and very interestingly, the same cortical areas (denoted by the same colors) were found to be present in the networks associated with 2 consecutive BNSs. This is typically the case for the left occipital pole (red color) and the left/right inferior occipital gyrus and sulcus (blue) that were identified at BNS1, BNS2, BNS3 and BNS4. Similarly, the left inferior temporal sulcus (cyan) remained present in networks associated with BNS3 and BNS4 while the anterior occipital sulcus (dark blue) was found to be active during BNS2 and BNS3.

These results suggest that the spatiotemporal brain dynamics of picture naming are characterized by transitions between overlapping networks associated with transiently stable states (from 30 ms to 160 ms), as summarized in Fig. 3B. Indeed, as depicted, four transitions (T1-T4) are characterized by a "continuous" flow of information processing as some network nodes common to consecutive BNSs were revealed by the proposed source connectivity approach. In addition, as detailed in the discussion, we could relate i) the timing of the 6 periods and ii) the identified networks associated with BNSs to well-established neuropsychological phases of the picture recognition/naming task (visual recognition, memory access and categorization, semantic processing, phonological encoding and articulation).

Discussion

Appropriate dense-EEG processing reveals the spatiotemporal dynamics of functional brain networks.

To our knowledge, this study is one of the very few attempts to characterize the spatiotemporal dynamics of brain networks over a short duration cognitive task second) from scalp EEG data. Results show that appropriate processing of dense-EEG recordings allows for identification of networks that are in agreement with brain regions involved in the same cognitive task and identified from other modalities (mainly fMRI and PET, [START_REF] Price | A review and synthesis of the first 20years of PET and fMRI studies of heard speech, spoken language and reading[END_REF], see below for details). However, these results go beyond those obtained with classical neuroimaging techniques, as the proposed method offers the unique advantage to track the network dynamics with high temporal (in the order of ms) and spatial (~1000 ROIs) resolution.

The good performance of this processing can be explained by three key steps to obtaining relevant networks, in terms of time and space features. The first step is the reconstruction of cortical sources distributed over a high-resolution mesh by solving the EEG inverse problem.

We previously realized a multifactor analysis to analyze the effect of the different factors that intervene in the EEG source connectivity analysis. This methodological study showed that the wMNE algorithm combined with dense electrode array (180 scalp electrodes) leads to optimal results [START_REF] Hassan | EEG Source Connectivity Analysis: From Dense Array Recordings to Brain Networks[END_REF]. The second step is the estimation of the functional connectivity using phase synchronization among gamma oscillations present in the time-course of reconstructed sources. Our results demonstrate that this step is crucial for identifying networks with high specificity in regard to the performed task. The third step is the segmentation, in time, of the cognitive process into brain network states (BNSs). Based on the k-means clustering of brain networks, we have developed an algorithm, originally at scalp EEG level at ms scale [START_REF] Mheich | A new algorithm for spatiotemporal analysis of brain functional connectivity[END_REF]. This algorithm was used here for the first time on cortical level networks. This segmentation procedure automatically leads to a timing that strongly matches the successive steps previously reported of brain processing from image perception to naming. This 'switching behavior' of the function connectivity networks has been very recently reported for resting states data using modeling approach [START_REF] Hansen | Functional connectivity dynamics: Modeling the switching behavior of the resting state[END_REF]. However, this time varying characteristics of the functional brain networks should be taken with caution as similar behaviors could be observed in random data (noise) as reported in Hlinka [START_REF] Hlinka | On the danger of detecting network states in white noise[END_REF], hence the importance of testing the results against null models which was performed in our paper.

Source functional connectivity vs. source localization

A major question that is addressed in this paper relates to the difference between the proposed network-based approach and the previously-used approach aimed at localizing the sources of activated regions during cognitive tasks [START_REF] Miozzo | Early parallel activation of semantics and phonology in picture naming: Evidence from a multiple linear regression MEG study[END_REF][START_REF] Salmelin | Dynamics of brain activation during picture naming[END_REF]. Although both methodologies (source-based and networkbased) lead to similar results regarding the involvement of the occipital and frontal lobes during the cognitive task, this study shows that the information extracted from dense-EEG recordings in both cases is dramatically different.

Conceptually, the fundamental difference between both approaches is that the source localization totally ignores all possible interactions between brain regions. When performing source localization analysis, the sources with highest amplitude (averaged at given time period or computed at the instant of peak amplitude of the signal) were classically retained [START_REF] Salmelin | Dynamics of brain activation during picture naming[END_REF]. However, to some extent (depending on threshold), this approach neglects the possible contribution of "low energy" sources. As shown in the supplementary figure S3, the threshold process can dramatically modify the results of the localization. For instance, when only the highest 50 or 75% amplitude sources are kept in the wMNE algorithm, the temporal lobe is not retained as being activated although it has been reported to participate in semantic processing in the exact same task [START_REF] Patterson | Where do you know what you know? The representation of semantic knowledge in the human brain[END_REF][START_REF] Price | A review and synthesis of the first 20years of PET and fMRI studies of heard speech, spoken language and reading[END_REF].

Conversely, the hypothesis behind the network-based approach is that sources can be synchronized regardless their amplitude. Our results show that edges represented in identified networks correspond to connections both among sources with high amplitude (occipital) and sources with low amplitude (temporal). We believe that the network-based approach allows revealing networks that are more specific to the performed task. An illustrative example is the absence of the dorsolateral prefrontal cortex in identified networks during the first 200 ms while this region is particularly active over the entire task in the source-localization approach, probably due to sustained attentional processes that are not strictly related to picture visualization, decoding and naming.

Brain networks involved in picture naming task

Regarding the successive steps of the picture naming task, our findings corroborate previous studies based on other modalities, mainly fMRI and PET (see [START_REF] Price | A review and synthesis of the first 20years of PET and fMRI studies of heard speech, spoken language and reading[END_REF] for review). For the first BNS (BNS1, 0:119 ms), results showed a network involving the inferior occipital, the lateral occipito-temporal sulcus and occipital pole. This period was shown to be related to the visual feature extraction preceding the object category recognition [START_REF] Thorpe | Speed of processing in the human visual system[END_REF][START_REF] Vanrullen | The time course of visual processing: from early perception to decisionmaking[END_REF]. Interestingly, the visual features obtained by Vanrullen et al. revealed response peaking at 120 ms after the onset which is very comparable to our limit for the first BNS (119 ms). For BNS2 (120:150 ms), the identified functional network mainly comprises the occipital regions (the bilateral inferior occipital, the left occipital pole, right anterior occipital and the left middle occipital and Lunatus). These regions are well known to play a capital role in the processing of visual information and object recognition [START_REF] Engel | Interhemispheric synchronization of oscillatory neuronal responses in cat visual cortex[END_REF][START_REF] Müller | Induced gamma-band responses in the human EEG are related to attentional information processing[END_REF]. Moreover, the gamma activity in this time period was shown to marker of object recognition and binding [START_REF] Martinovic | Induced gamma-band activity elicited by visual representation of unattended objects[END_REF][START_REF] Müller | Induced gamma-band responses in the human EEG are related to attentional information processing[END_REF].

For BNS3 (151:190 ms), results also indicate a mainly occipital network but with an implication of the bilateral inferior temporal sulcus. This system is known to be related to lexical retrieval, lemma retrieval and lemma selection [START_REF] Indefrey | The spatial and temporal signatures of word production components[END_REF]. It is also involved in semantic working memory system when someone tries to remind the name of the objects [START_REF] Martin | Semantic memory and the brain: structure and processes[END_REF]. In their study, the authors show a discrepancy in the temporal lobe involvement for objects versus animals with more activity in the inferior temporal sulcus for objects and in the superior temporal sulcus for animals. Our picture set comprises 39 animals versus 109 objects or non-animal images. The dominant representation of objects in our experimental set could have shaped this part of the graph favoring the inferior temporal sulcus at the expense of the superior temporal sulcus. During BNS4 (191:320 ms), the network involves the left inferior temporal gyrus in addition to the inferior temporal sulcus.

These regions were stated to be in direct relation to semantic processing [START_REF] Martin | Semantic memory and the brain: structure and processes[END_REF]. It is also the time window in which the N200 classically appear. The N200 is a marker of semantic processing in go/no-go tasks [START_REF] Thorpe | Speed of processing in the human visual system[END_REF]. Together with the appearance at this stage of frontal nodes, we assume that this large BNS could also integrates the access to phonological forms during overt naming which has been shown by [START_REF] Graves | A neural signature of phonological access: distinguishing the effects of word frequency from familiarity and length in overt picture naming[END_REF]. In a somewhat interesting way, we found that this BNS ends 10 ms before the syllabification step in the model of Indefrey and Levelt [START_REF] Indefrey | The spatial and temporal signatures of word production components[END_REF].

At BNS5 (321:480 ms), the functional network involves the left superior precentral and the right postcentral sulcus along with the left orbital sulci and the left superior insula. This network, near to the sensory-motor cortex was reported to be engaged into phonetic and articulatory process of speech [START_REF] Dronkers | A new brain region for coordinating speech articulation[END_REF][START_REF] Levelt | An MEG study of picture naming[END_REF]. Finally, for the last BNS (BNS6, 481:535 ms), the network was found at the left insular gyrus, the left inferior insular circular sulcus and the right orbital sulcus. This network is typically associated with the beginning of the naming process and the speech articulation as well as listening to own speech (external self-monitoring) [START_REF] Indefrey | The spatial and temporal signatures of word production components[END_REF][START_REF] Levelt | An MEG study of picture naming[END_REF]. At this stage, the variability between subject's response latencies has become incredibly strong and those results should be considered with caution. A backward analysis of the EEG records and the corresponding BNSs (not presented here) can give us more precise results about the last graphs.

The occipital cortex remains active during several BNSs.

Strikingly, our results indicated that the occipital visual areas (especially the inferior occipital area) remain active for 400 ms, over the first four BNSs. This is likely not an artifact due to open eyes for instance because left occipital region was found to be inactive for the last BNS whereas the pictures are still present on the screen. Also, it cannot be due to a variable delay in the behavioral response as the timing of the very early steps of the occipital activity for object identification is only very slightly modified due to attentional effects and the N1 is not modified at all [START_REF] Mangun | Sustained visual-spatial attention produces costs and benefits in response time and evoked neural activity[END_REF]. Those regions are part of what the literature calls the lateral occipital complex (LOC) [START_REF] Grill-Spector | The lateral occipital complex and its role in object recognition[END_REF]. They are particularly involved in object representation and for analyzing objects shape independently of other low level characteristics such as color texture. Those characteristics not only serve to address the mental lexicon but also form a part of the large-scale network that builds the semantic representation of objects [START_REF] Grill-Spector | The lateral occipital complex and its role in object recognition[END_REF][START_REF] Huth | A continuous semantic space describes the representation of thousands of object and action categories across the human brain[END_REF]. Those representations are held in visual working memory during the complete process [START_REF] Harrison | Decoding reveals the contents of visual working memory in early visual areas[END_REF]. It could then take part of a large-scale network characterized by high gamma synchronization contributing to information processing during almost the whole cognitive process of naming.

This persistent occipital cortex activation suggests that the ventral stream does not simply correspond to a sequential Feed-Forward activation [START_REF] Gilbert | Top-down influences on visual processing[END_REF] of involved brain regions but instead, to a resonance process within a large-scale network. From these results, we hypothesize that, in this network, the image is first processed in the occipital areas of the visual cortex, inducing a sustained gamma activity that lasts for on average 400 ms. This oscillatory activity is then progressively and sequentially compared to oscillations in the other areas. Finally, image categorization, memory access and semantic processing emerge from gamma synchronization among involved brain regions. These results are also in accordance with the resonance properties of specific networks of neurons that are able to selectively respond to inputs at preferred frequencies [START_REF] Akam | Oscillatory multiplexing of population codes for selective communication in the mammalian brain[END_REF][START_REF] Hutcheon | Resonance, oscillation and the intrinsic frequency preferences of neurons[END_REF].

Methodological considerations

In the present study, the whole processing was performed on group-averaged data. The connectivity matrix was computed using the phase locking value (PLV) method at each time instant for each subject and all matrices were averaged over the 21 subjects, yielding a group-averaged connectivity matrix at each time sample of the post-stimulus response. The segmentation algorithm was applied to this averaged data. The main advantage of this approach is that it preserves the networks common to all subjects of the group and reduces the inter-individual variability. Consequently, results obtained from this analysis represent the spatiotemporal dynamics of networks that occur most constantly within the group of subjects.

It is also noteworthy that we assumed i) that there is an anatomical correspondence among subjects (a template 3D mesh was used as the source model) and ii) that the transitions between the functional connectivity states occur in a repeatable manner over subjects.

However, the inter-subject variability was not totally ignored in our analysis. The index "network presence" (see supplementary materials) of each BNS was computed among all the subjects. The results show that the first three BNSs have the highest network presence (80%, 82% and 81% for BNS1, BNS2 and BNS3 respectively). Conversely, the network presence values decreased after BNS3 to reach 64%, 62% and 41% for BNS4, BNS5 and BNS6 respectively. We are aware that due to strong behavioral reaction time discrepancies, the resonance duration of the nodes within the occipital areas over the BNS4 and BNS5 depends on the variability of the cognitive steps that follows the object identification. Those cognitive steps are importantly modulated by attentional resources and the 400 ms resonance duration before the BNS6 could be due to a mix of long BNS4 in some individuals or for some specific events with short BNS5 of other faster subjects.

Regarding the functional connectivity measure, we used the PLV method [START_REF] Lachaux | Measuring phase synchrony in brain signals[END_REF] which provides high performance to detect the inter-trial synchronization at each time instant. However, the PLV method requires a relatively high number of stimuli (148 in our case for each subject) to be correctly estimated. Consequently, this method cannot be easily applied to on-going activity, as for the analysis of resting state networks for instance.

Although EEG source connectivity reduces the problem of field spread to a certain extent, it

does not yet, provide a perfect solution. The field spread effect is an open question and unsolved issue and it is indeed one of the main challenges when performing EEG/MEG inverse solution. In the connectivity context, the main effect of the field spread is referred to 'source leakage', which denotes the possible 'artificial' correlations that can exist between very close sources. In our study, we have considerably limited this problem. We firstly decided to use PLV, which has been shown to be very efficient to reduce field spread (REF).

The second step is the 'clustering' or the anatomical parcellation, which consist of reconstructing the sources on a high spatial resolution cortex mesh then averaging the dynamics sources located in the same ROIs. The functional connectivity was then computed between these averaged sources. The averaging over ROIs increases the distance between the positions of the sub-regions, which reduces some artificial correlations between very close sources (at the same ROIs). In addition, our results showed high number of long range connections. Very few approaches have been proposed recently to deal with the source leakage by either normalizing the edges weights by the distance between the nodes or removing the very close edges. However, each of these approaches has its advantage and disadvantages and can remove 'real' connections in most cases (REF Shofflen 2009).

Nevertheless, it is important to keep in mind that the source reconstruction algorithms can reduce the volume conduction problem but not address it completely.

Conclusion

We used dense-EEG recordings during a picture naming task to characterize the spatiotemporal dynamics of functional brain networks. We showed that the picture naming task can be divided into six brain network states (BNSs) characterized by significantly high gamma synchronization. Results revealed that fast transitions occur between these BNSs and last from 30 ms to 160 ms. The spatial location of nodes and edges in the identified networks, as well as the precise timing of transitions, show a precise cognitive function (visual recognition, semantic processing and speech processing) can be associated with each network.

In addition, networks associated with BNSs partially overlap. We hypothesized that the persistence of several brain regions in successive BNSs participates to fast and efficient information in the brain.

Figure 1:

Structure of the investigation: A) Dense EEG (256 electrodes) were recorded during picture naming task. Structural MRI images were segmented and anatomically parcellated [START_REF] Destrieux | Automatic parcellation of human cortical gyri and sulci using standard anatomical nomenclature[END_REF] using Freesurfer (Fischl, 2012) to obtain 148 regions. These 148 regions were subdivided using Brainstorm [START_REF] Tadel | Brainstorm: a user-friendly application for MEG/EEG analysis[END_REF] to obtain higher spatial resolution (about 1000 regions of interests). The inverse problem was then solved using the weighted Minimum Norm Estimate algorithm. The time series of the reconstructed sources were obtained. The functional connectivity between the reconstructed sources was computed using the Phase Locking Value method. A high resolution functional connectivity matrix was obtained and the corresponding functional brain network was visualized. B) This procedure was carried out at millisecond scale and a high number of functional connectivity matrices were obtained. A Kmeans clustering approach was used to obtain the brain network states (BNSs) using a segmentation algorithm [START_REF] Mheich | A new algorithm for spatiotemporal analysis of brain functional connectivity[END_REF] summarized in materials and methods section and detailed in the supplementary document.

  According to the linear discrete equivalent current dipole model, EEG signals S(t) measured from Q channels can be expressed as linear combinations of P time-varying current dipole sources D(t): S= G. D+B

  xy t represents the general term of the average adjacency matrix PLV( ) t which defines a functional connectivity network N at each time t,

Figure 3 :

 3 Figure 3:

Figure 4 :

 4 Figure 4:

Figure 5 :

 5 Figure 5:

  used MATLAB ® (2007a, MathWorks Inc.) in the entire process: from EEG preprocessing, source reconstruction, functional connectivity analysis, and computation of the graph parameters and visualization of the brain networks. These processing steps were carried out by using several tools mainly Freesurfer http://surfer.nmr.mgh.harvard.edu/ for surface parcellation, Brainstorm http://neuroimage.usc.edu/brainstorm/ for source reconstruction, EEGLAB http://sccn.ucsd.edu/eeglab/ for EEG preprocessing, and other homemade codes.3D brain networks were visualized using a recently developed tool in our team called EEGNET. The network measures were computed from the brain connectivity toolbox (BCT)

https://sites.google.com/site/bctnet/ that was integrated in EEGNET.

Acknowledgements

This work has received a French government support granted to the CominLabs excellence laboratory and managed by the National Research Agency in the "Investing for the Future" program under reference ANR-10-LABX-07-01. It was also financed by the Rennes University Hospital (COREC Project named conneXion, 2012-14). This work was also supported by the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement n° 290901. The authors also want to thanks Deok-Hee Kim-Dufor for her help in analyzing the voice onset time of the recorded oral responses.