Anis Amziane
email: anis.amziane@univ-lille.fr

Olivier Losson
email: olivier.losson@univ-lille.fr

Benjamin Mathon
email: benjamin.mathon@univ-lille.fr

Aurelien Dumenil
email: a.dumenil@somme.chambagri.fr

Ludovic Macaire
email: ludovic.macaire@univ-lille.fr

Aurélien Dumenil

Frame-based reflectance estimation from multispectral images for weed identification in varying illumination conditions

Keywords: Multispectral imaging, Reflectance estimation, Illumination, Precision farming, Weed detection, Linescan camera

est destinée au dépôt

I. INTRODUCTION

In most cases, removal of the weeds in agricultural fields involves the application of abundant quantities of chemical herbicides, which is harmful to the environment regardless of their success in increasing crop productivity. To optimize the application of herbicides in crop fields, precision spraying can be considered thanks to recent advances in imaging sensors. The sensors originally used in agricultural fields are either monochromatic or color-based ones [START_REF] Mulla | Twenty five years of remote sensing in precision agriculture: Key advances and remaining knowledge gaps[END_REF]. A critical limitation is that they use average spectral information over wide wavelength ranges, resulting in a lack of detailed information in specific narrow bands. During the last decade, more sophisticated multispectral sensors have been manufactured and deployed in crop fields leading to precision farming applications like pest and disease detection [START_REF] Mahlein | Detection, identification, and quantification of fungal diseases of sugar beet leaves using imaging and non-imaging hyperspectral techniques[END_REF], plant classification [START_REF] Haug | Plant classification system for crop /weed discrimination without segmentation[END_REF], and weed detection [START_REF] Wendel | Self-supervised weed detection in vegetable crops using ground based hyperspectral imaging[END_REF]- [START_REF] Lin | Detection of corn and weed species by the combination of spectral, shape and textural features[END_REF]. Multispectral cameras collect spectral data over a wide spectral range and provide the ability to investigate spectral responses of soils and vegetated surfaces in narrow spectral bands. Two categories of devices can be distinguished in multispectral image acquisition [START_REF] Hagen | Review of snapshot spectral imaging technologies[END_REF]. "Snapshot" (multi-sensor or filter arraybased) devices build the image from a single shot. Although this technology provides multispectral images at video frame rate, the few acquired channels may not be sufficient to fully explore the vegetation spectral signatures. "Multishot" (tunable filter or illumination-based, push-broom, spatio-spectral linescan) devices build the image from several and successive frame acquisitions. Although restricted to still scenes, they provide images with a high spectral resolution. In this study, we use a multishot camera to acquire outdoor multispectral radiance images of plant parcels in a greenhouse. From this radiance information, the spectral reflectance is estimated as an illumination-invariant spectral signature of each species.

Several methods have been proposed to compute reflectance thanks to prior knowledge about cameras or illumination conditions [START_REF] Shen | Reflectance reconstruction for multispectral imaging by adaptive Wiener estimation[END_REF]- [START_REF] Heikkinen | Evaluation and unification of some methods for estimating reflectance spectra from RGB images[END_REF]. In field conditions, typical methods first estimate the illumination by including a reference device (a white diffuser or a ColorChecker Chart) in the scene [START_REF] Feyaerts | Multi-spectral vision system for weed detection[END_REF], [START_REF] Bourgeon | Field radiometric calibration of a multispectral on-the-go sensor dedicated to the characterization of vineyard foliage[END_REF]- [START_REF] Uto | Characterization of rice paddies by a UAV-mounted miniature hyperspectral sensor system[END_REF]. Then, the reflectance is estimated at each pixel p by channel-wise dividing the value of the radiance image at p by the pixel values characterizing the white diffuser or the ColorChecker white patch. In [START_REF] Khan | Illuminant estimation in multispectral imaging[END_REF], an extension to the multispectral domain of four algorithms traditionally applied to RGB images is proposed to estimate the illumination. In [START_REF] Eckhard | Outdoor scene reflectance measurements using a bragg-grating-based hyperspectral imager[END_REF], a bragg-grating-based multispectral camera is used to acquire outdoor radiance images. Two white diffusers are used for reflectance computation, one is included at the bottom image border and the other is fully visible. In [START_REF] Alexander | Illumination compensation in ground based hyperspectral imaging[END_REF], illumination variation is compensated in outdoor condition using a subspace model proposed in [START_REF] Mark | Analytic solution for separating spectra into illumination and surface reflectance components[END_REF], and enhanced by learning a historical database of illumination and reflectance spectra along various natural lighting conditions. A multispectral camera is used in conjunction with a skyward pointing spectrometer in [START_REF] Chuiqing | Fusion of multispectral imagery and spectrometer data in UAV remote sensing[END_REF] to estimate reflectance from the acquired scene radiance. These methods require specific devices and knowledge about the spectral sensitivity functions (SSFs) of the sensor filters, and may include a learning step. They also often assume constant incident illumination during the successive frame acquisitions (scans). In outdoor conditions however, the illumination may vary significantly across scans that last several seconds (e.g., 12s with our camera for an integration time of 1.5ms Irradiance (W

• m -2 • nm -1 )
Acquisition time (s)

Fig. 1: Illumination variation at wavelength 500nm during an acquisition of an outdoor multispectral image.

that is suited to most outdoor conditions). Fig. 1 illustrates this phenomenon by showing the variation of the incident irradiance at wavelength 500nm during one of our outdoor frame acquisitions, measured by the AvaSpec-ULS2048L spectrometer. In this paper, we propose a learning-free reflectance estimation method that is robust to illumination variations during the acquisition of frames used to provide multispectral images.

The paper is organized as follows. In Sec. II, we give details about the multispectral image acquisition process in outdoor conditions using a linescan camera. In Sec. III, we propose a method for scene reflectance estimation from a multispectral radiance image acquired in uncontrolled and varying illumination conditions. Sec. IV presents an experimental evaluation of the proposed method, and weed/crop segmentation results using the estimated reflectance. Finally, conclusions are drawn in Sec. V.

II. MULTISPECTRAL IMAGE ACQUISITION

A. Image formation model

We consider a multispectral camera that embeds a sensor covered with spectrally sensitive filters and provides a multispectral radiance image I = {I k } K k=1 with K spectral channels thanks to the acquisition of frames. Assuming homogeneous spectral sensitivity of the sensor, the value I k p of channel I k at pixel p can be expressed as:

I k p = Q τ • Ω E(λ) • R p (λ) • T k (λ) • o p (λ) dλ , ( 1
)
where τ is the frame integration time, Ω is the working spectral domain, E(λ) ∈ [0, 1] is the relative spectral power distribution of the illumination which is assumed to homogeneously illuminate all surface elements of the scene, R p (λ) ∈ [0, 1] is the spectral reflectance at wavelength λ of the surface element s observed by pixel p, o p (λ) ∈ [0, 1] is the optical attenuation of the lens for pixel p at λ. The radiance

E(λ) • R p (λ) • o p (λ)
is filtered according to the SSF T k (λ) of the filter k associated to a narrow spectral band. The value I k p is given by the quantization of the received energy on B bits by the function Q.

B. Snapscan camera

The Snapscan [START_REF] Pichette | Fast and compact internal scanning CMOS-based hyperspectral camera: the Snapscan[END_REF] is a multispectral camera manufactured by IMEC that embeds a single matrix sensor, covered by a series of narrow stripes of Fabry-Perot integrated filters. It provides K = 141 spectral channels whose central wavelengths range from λ 1 = 475.1nm to λ 141 = 901.7nm with a variable center step (from 0.5nm to 5nm). Specifically, each filter is associated to 5 adjacent rows of 2048 pixels, and samples a band from the visible (VIS) or near infrared (NIR) spectral domain according to its SSF T k (λ) with a full width at half maximum between 2nm and 10nm. During acquisitions, both the object and camera remain static. The sensor moves linearly from bottom to top behind the lens inside the Snapscan with a spatial step of 5 pixel rows, and records a two-dimensional frame f i at each time t i , i = 1 . . . N . A frame f i is spatially organized as juxtaposed stripes of 5 consecutive pixel rows. Each stripe f k i , k = 1 . . . K, contains the spectral information of the scene radiance filtered according to the SSF T k (λ) of the filter k centered at wavelength λ k . The multispectral image should then be assembled from N successively acquired frames {f i } N i=1 (see Fig. 2). All the stripes associated to filter k in all the acquired frames are stacked to provide an assembly {f k i } N i=1 . The scene part common to all the assemblies is considered to form the spectral channels. Finally, the Snapscan delivers a multispectral radiance image I = {I k } K k=1 composed of K = 141 channels quantized on B = 10 bits. Note that before the frame acquisitions, the Snapscan uses its internal shutter to acquire a dark frame f dark whose values are subtracted pixel-wise from the acquired frames. Therefore, the pixel value expressed by ( 1) is assumed to be free from thermal noise.

III. REFLECTANCE ESTIMATION

A. Reflectance computation

To compute reflectance according to the linear image formation model of (1), the following assumptions need to be verified: (i) The illumination is spatially uniform and does not vary during the successive acquisition of the frames. (ii) Each of the Fabry-Perot filters has an ideal SSF:

T k (λ) = δ(λ -λ k ) = 1 if λ = λ k , 0 otherwise. ( 2
)
(iii) Each SSF of the Fabry-Perot filters is angle-independent, thus vignetting effect only depends on pixel location. The discrete version of (1) can be expressed as:

I k p = τ • K l=1 E(λ l ) • R p (λ l ) • T k (λ l ) • o p (λ l ). (3)
Note that the quantization function Q is omitted here since the different terms are already considered as quantized.

t 1 t 2 t 3 t 4 t N f 1 1 f 1 1 f K 1 f 1 2 f 1 2 f K 2 f 1 3 f 1 3 f K 3 f 1 4 f 1 4 f K 4 f 1 N f 1 N f K N Acquisition time Illumination cond.
Scanning direction 2048 pixels 2048 pixels 5 pixels

Snapscan sensor

Scanned object

Stripe stack (first filter)

5N pixels Assumption (ii) is not verified in practice: there is no multispectral camera with such a fine resolving power. However, to the extent of our knowledge, the filters embedded in the Snapscan are very narrow and approximate this assumption to the best at present. We can then write:

I 1 p I K p {f 1 i } N i=1 assembly Discarded part of {f 1 i } N i=1 assembly Discarded part of {f K i } N i=1 assembly I 1 I K Wavelength (nm) λ 1 λ K x Assembled image I
I k p = τ • E(λ k ) • R p (λ k ) • o p (λ k ). (4)
The reflectance of a surface element s observed by pixel p for the spectral band centered at λ k is then expressed as:

R p (λ k ) = I k p τ • E(λ k ) • o p (λ k ) . (5)

B. Reflectance estimation under constant illumination

The term E(λ k )•o p (λ k ) in ( 5) characterizes the illumination in the sensor domain for channel k, and can be experimentally estimated by acquiring a full-field image I wd of a white diffuser. This calibration tile is supposed to reflect the incident light for all λ ∈ Ω with a constant diffuse reflection factor ρ wd (ρ wd = 95% for our Sphere Optics Spectralon ® ). Using the same image formation model and assumptions as for I (see (4)), each value of I wd is then expressed as:

I k wd,p = τ wd • E(λ k ) • ρ wd • o p (λ k ), (6)
where τ wd is the frame integration time used to acquire I wd . Plugging (6) into (5) then yields the estimated reflectance Rk p for channel k at pixel p:

Rp (λ k ) def = Rk p = ρ wd • I k p I k wd,p • τ wd τ . (7)
Note that this reflectance estimation method implicitly compensates the vignetting effect. Since the white diffuser and the object occupy the same (full) field of view, I k p and I k wd,p are affected by the same optical attenuation whose effect vanishes after division. We consider the reflectance estimated by this method as a reference and denote it as Rref .

C. Reflectance estimation under varying illumination

The Snapscan assembles a multispectral radiance image from individual frames acquired at different successive times. Due to illumination variations during frame acquisitions in outdoor (see Fig. 1), the illumination conditions associated to the values of the assembled radiance image vary both across pixel rows for a given channel and across channels for a given row (see Fig. 2). To estimate the reflectance, we then need to estimate the illumination condition for each channel and each row.

Inspired by [START_REF] Eckhard | Outdoor scene reflectance measurements using a bragg-grating-based hyperspectral imager[END_REF], [START_REF] Alexander | Illumination compensation in ground based hyperspectral imaging[END_REF], we propose to mount a white diffuser on the acquisition device so that the sensor vertically observes a portion of it at the right border of each frame. Once the N frames have been assembled as an image, an area W of width 250 pixels represents the white diffuser (see Fig. 3). Each row of this area is used to estimate the specific illumination condition under which the scene part associated to the same pixel row has been acquired.

Because reflectance cannot be estimated according to (7) any longer, the vignetting effect cannot be compensated as in Sec. III-A. To accurately estimate the reflectance, we need to spatially correct the vignetting effect on each channel of the radiance image I.

1) Vignetting correction: To compute a correction factor for each pixel, we use a look-up table (LUT) approach [START_REF] Yu | Practical anti-vignetting methods for digital cameras[END_REF] because it requires no knowledge about the optical device behavior. To generate LUT correction factors, one first acquires a full-field image I wd of a white diffuser under a spatially uniform and constant illumination, and in a dark room to avoid ambient light influence. The vignetting effect in the object radiance image I is then corrected pixel-wise and channelwise to provide the corrected image Ĩ as:

Ĩk p = I k wd I k wd,p • I k p . ( 8
)
The first term is the correction factor. It implies I k wd , namely the median value of m pixels with the highest values over I k wd to discard saturated or defective pixel values.

2) Row-wise (rw) based reflectance estimation: Since pixels in each row of Ĩ are associated to the same illumination condition according to the spatial uniformity assumption, we estimate the reflectance from the radiance image Ĩ in a rowwise manner (see Fig. 3). Moreover, the attenuation is spatially uniform (and equal to a constant α k for all p associated to channel index k) after vignetting correction. The reflectance at pixel p associated to channel index k is then deduced from (5) as:

R p (λ k ) = Ĩk p τ • E y(p) (λ k ) • α k . ( 9
)
For a given row y(p), the illumination E y(p) (λ k ) is the same at any pixel p. We estimate it from the median value Ĩk

W,y(p)
of the m highest pixel values of the white diffuser (W) in this row. This median value satisfies [START_REF] Lin | Detection of corn and weed species by the combination of spectral, shape and textural features[END_REF], which provides the rw-based estimated reflectance as:

Rk rw,p = ρ wd • Ĩk p Ĩk W,y(p) . ( 10
)

IV. EVALUATION OF REFLECTANCE ESTIMATION A. Overview

An image acquisition campaign was conducted in a greenhouse under skylight by the Chambre d'Agriculture de la Somme in early April 2019. The targeted plants are beet, wheat, bean, and barley for crops, and datura, thistle, goosefoot, foxtail, and barnyard grass for weeds. For our exper- We compare the quality of the rw method (setting m = 11) with the results obtained by two methods that are generally dedicated to images acquired with snapshot or linescan cameras in controlled illumination conditions. The classical whiteaverage (wa) method uses the surface of a white diffuser included in the scene to estimate the illumination for each channel. We adopt the wa implementation by considering a subset W s of W in our images (see dashed green rectangle in Fig. 3). The max-spectral (ms) method provides good performance when estimating the illumination in multispectral images with various band combinations [START_REF] Khan | Illuminant estimation in multispectral imaging[END_REF]. It uses the maximum pixel value within each channel. The reflectance is estimated at each pixel in each channel by wa and ms methods as:

Rk wa,p = ρ wd • Ĩk p 1 |W s | p∈Ws Ĩk p , (11)
and:

Rk ms,p = Ĩk p max p∈X Ĩk p , (12)
where |W s | is the number of white diffuser pixels considered for wa implementation, and X is the set of vegetation and soil pixels (i.e., all image pixels but W and those of the ColorChecker).

B. Reflectance estimation error

To evaluate the accuracy of reflectance estimation, we use the patches of the ColorChecker (see Fig. 3). Let R * , * ∈ {rw, wa, ms}, denote the reflectance image estimated by either of the three tested methods (see ( 10)-( 12)). All surface elements of a patch are assumed to have the same spectral response. Each patch P j , j = 1 . . . 24, is then represented by a K-dimensional reflectance vector R * ,Pj whose k th component is computed as the average value in channel k over all the pixels characterizing the patch. This vector is compared to the reference reflectance Rref,Pj of the same patch estimated according to [START_REF] Hagen | Review of snapshot spectral imaging technologies[END_REF] from an image of the ColorChecker acquired in laboratory under controlled illumination conditions. Note that the ColorChecker is placed at various spatial coordinates in the images of our database. In case of high vegetation density, it is placed on the top of a wooden block to prevent patch occlusions by leaves. Thus we use a 15 × 15 pixel window around the automatically detected center of each patch to ensure all these pixels belong to the same patch whatever the image.

We objectively asses each estimated reflectance image thanks to the mean absolute error (MAE) and the angular error ∆θ of each patch P j , given by:

MAE( Rref,Pj , R * ,Pj ) = 1 K K k=1 | Rk ref,Pj -Rk * ,Pj |, (13)
and: ∆θ( Rref,Pj , R * ,Pj ) = arccos Rref,Pj , R * ,Pj Rref,Pj

2 • R * ,Pj 2 , ( 14
)
where , is the dot vector product, and • 2 is the Euclidean norm. When ∆θ of two vectors (spectra in our case) is equal to zero, it means that these two vectors are collinear. We compute the average mean absolute error MAE * and angular error ∆θ * over all patches of all estimated reflectance images. We also compute the average channel-wise mean absolute error MAE k * over all patches of all reflectance images.

Tab. I shows the average mean absolute and angular errors over all images in the database for the three tested methods. We can see that in terms of quality of reflectance estimation, the rw-based method provides the lowest average mean absolute and angular errors. Considering the channel-wise error (see Fig. 4), the highest error rates in the range [475nm, 750nm] (VIS domain) are obtained by ms because it only analyzes pixels of background and vegetation that strongly absorb the incident light. Hence, illumination estimation is biased using ms in this domain. In the range [750nm, 901nm] (NIR domain), the light reflection factor of vegetation is high (it appears brighter in Fig. 3), and the ms method gives very similar results as rw. The wa method gives acceptable results along the working domain [475nm, 901nm] of the Snapscan camera, while our rw method provides the lowest error rates.

C. Multispectral image segmentation

We now briefly evaluate the contribution of the rw-based reflectance estimation method for segmentation tasks. To assess the robustness of rw against illumination variations, we consider an extra outdoor radiance image of a crop (beet) and weed (thistle) species (see Fig. 5), acquired on a different day than those present in our database. Vegetation pixels are separated from the background using the NDVI index [START_REF] Thenkabail | Evaluation of narrowband and broadband vegetation indices for determining optimal hyperspectral wavebands for agricultural crop characterization[END_REF]. A pixel p is considered as a vegetation pixel if its NVDI index is greater than a manual threshold T , namely: R140 rw,p -R68

V. CONCLUSION

This paper proposes a method to estimate spectral reflectance from multispectral radiance images acquired in outdoor illumination conditions. We show that classical methods of reflectance estimation are not adapted in these conditions for multishot devices such as the Snapscan. Indeed, the illumination that varies during the frame acquisitions affects the reflectance estimation quality. Therefore, we propose to estimate the illumination at a frame level following the proposed row-wise (rw) approach. As an application of the rw-based reflectance, an image segmentation by pixel classification is performed. As a future work, the segmentation process will be extended to more images with complex scenes. Furthermore, the most suited spectral bands for this task should be identified and used to specify a multispectral camera that is less sensitive to outdoor illumination variation.

Fig. 2 :

 2 Fig. 2: Frame acquisition and assembly in varying conditions (top) and K-channel multispectral image (bottom).

Fig. 3 :

 3 Fig. 3: rw-based reflectance estimation (k = λ 100 = 780nm).

Fig. 4 :

 4 Fig. 4: Channel-wise MAE of reflectance estimation.

Fig. 5

 5 show the segmentation results obtained on the input multispectral image rendered as RGB in Fig.5(a) using the CIE Standard Illuminant D65 with the CIE 1931 2 • Standard Observer. Segmentation using radiance spectra (see Fig. 5(c)) leads to thistle over-detection whereas waand msbased reflectance spectra (see Figs. 5(e) and 5(f)) lead to thistle under-detection. rw-based spectra (see Fig. 5(d)) provide a good discrimination between the two classes.

5 :

 5 Vegetation segmentation results. Beet is shown as green and thistle as red.

TABLE I :

 I REFLECTANCE ESTIMATION ERRORS

	Method * MAE * (%) ∆θ * (rad)
	rw		7.479	0.092
	wa	10.188	0.094
	ms	20.636	0.354
	40
				rw
	20 30 * (%) MAE k			wa ms
	10
		VIS domain		NIR domain
	0
	500	600	700	800	900
		Wavelength (nm)

TABLE II :

 II NUMBER OF PIXELS PER CLASS

	Class Extracted pixels Balanced set Test pixels
	Beet	26,000	20,000	674,131
	Thistle	20,000	20,000	340,866
	Total	46,000	40,000	1,014,997
				rw,p rw,p + R68 R140 rw,p	≥ T ,	(15)
				where λ 68 = 678.2nm and λ 140 = 899.2nm with the
				Snapscan. T is set to 0.45 in our case. Note that a vegetation pixel in any rw-based reflectance image is also considered
				as such in its corresponding radiance image and reflectance
				images estimated by wa and ms.
				The vegetation pixels are manually annotated by an expert
				in agronomy that builds the segmentation ground truth as two
				classes. Our radiance database contains 13 beet images and 10
				thistle images with different growth stages. From each image
				we randomly extract 2000 pixels, each being characterized
				by a feature vector of reflectance (or radiance) values of
				dimension K = 141. To reduce residual noise and within-
				class variability, the reflectance/radiance is averaged channel-
				wise over a 5 × 5 pixel window around each pixel. To analyze equiprobable classes, the pixel sets are balanced against the
				lowest class size (see Tab. II).

The segmentation process obeys a supervised scheme. All pixels in the balanced set are used as learning samples. We use the non-parametric 5-nearest neighbors classifier and the Euclidean distance measure.

Tab. III shows the classification accuracies. The rw method

TABLE III :

 III CLASSIFICATION ACCURACIES .83 62.91 provides the highest classification accuracy, which suggests that rw-based reflectance is the most robust against illumination variations.

	Spectral signature radiance	rw	wa	ms
	Accuracy (%)	76.86	88.93 77

ACKNOWLEDGEMENT

This work was supported by the Région Hauts-de-France, the Chambre d'Agriculture de la Somme, and the ANR-11-EQPX-23 IrDIVE platform.


Abderezak Salmi
email: salmi-abderezak@hotmail.fr

Kamal Hammouche
email: kamal_hammouche@yahoo.fr

Ludovic Macaire
email: ludovic.macaire@univ-lille.fr

Similarity-Based Constraint Score For Feature Selection

Keywords: Constraint score, Feature selection, Pairwise constraints, Similarity matrix

come

Introduction

In machine learning and pattern recognition applications, such as data mining and image analysis, datasets are often characterized by a large number of features. The processing of such high-dimensional data requires large memory storage and high computational time, and may lead to poor learning performance [START_REF] Kalakech | Constraint scores for semisupervised feature selection: A comparative study[END_REF][START_REF] Sheikhpour | A survey on semi-supervised feature selection methods[END_REF]. To address these drawbacks, the dimensionality of data is often reduced by using feature selection to remove redundant or irrelevant features.

Typically, feature selection methods can be categorized into three types: filter, wrapper, and embedded methods [START_REF] Sheikhpour | A survey on semi-supervised feature selection methods[END_REF][START_REF] Cai | Feature selection in machine learning: A new perspective[END_REF]. Filter methods evaluate features independently of the classification algorithm, while wrapper methods exploit a classification algorithm to evaluate the relevance of features. Embedded methods embed feature selection into the learning algorithm. Because filter methods are not dependent on any classification scheme, they have better generalization ability than wrapper and embedded methods [START_REF] Sheikhpour | A survey on semi-supervised feature selection methods[END_REF][START_REF] Cai | Feature selection in machine learning: A new perspective[END_REF].

According to the availability of prototypes (i.e., labeled data samples that represent classes), feature selection methods can also be divided into unsupervised, supervised, and semi-supervised approaches [START_REF] Kalakech | Constraint scores for semisupervised feature selection: A comparative study[END_REF][START_REF] Sheikhpour | A survey on semi-supervised feature selection methods[END_REF][START_REF] Cai | Feature selection in machine learning: A new perspective[END_REF]. Supervised feature selection only uses prototypes to measure the correlation of each feature with the class labels, while unsupervised feature selection uses only unlabeled data samples to evaluate the feature capacity to preserve the intrinsic data structure [START_REF] Kalakech | Constraint scores for semisupervised feature selection: A comparative study[END_REF]. Semi-supervised feature selection takes into account both prototypes and unlabeled data samples to evaluate the relevance of features.

In supervised and semi-supervised learning frameworks, besides class labels of prototypes, the available information can be also expressed by must-link and cannot-link constraints. A must-link constraint specifies that two data samples belong to the same class, while a cannot-link constraint specifies that two data samples belong to different classes [START_REF] Wagstaff | Constrained k-means clustering with background knowledge[END_REF]. Pairwise constraints can be provided by the user or easily generated from a small number of prototypes.

Must-link and cannot-link constraints are used to estimate the relevance of features via score functions, called constraint scores [START_REF] Kalakech | Constraint scores for semisupervised feature selection: A comparative study[END_REF][START_REF] Sheikhpour | A survey on semi-supervised feature selection methods[END_REF]. Zhang et al. [START_REF] Zhang | Constraint score: A new filter method for feature selection with pairwise constraints[END_REF] proposed two supervised constraint scores that use only pairwise constraints to evaluate the relevance of features. Zhao et al. [START_REF] Zhao | Locality sensitive semi-supervised feature selection[END_REF] defined a semi-supervised constraint score that analyzes both pairwise constraints and unlabeled data samples for feature selection. Kalakech et al. [START_REF] Kalakech | Constraint scores for semisupervised feature selection: A comparative study[END_REF] combined an unsupervised score com-puted from unlabeled data samples with a supervised score that is computed from the pairwise constraints. This score is predicted to be less sensitive to constraint changes. Two semi-supervised constraint scores that assess the ability of a feature to preserve the local properties of unlabeled data samples while respecting pairwise constraints, have been proposed by Benabdeslem et al. in [START_REF] Benabdeslem | Constrained laplacian score for semisupervised feature selection[END_REF] for the former and in [START_REF] Benabdeslem | Efficient semi-supervised feature selection: constraint, relevance, and redundancy[END_REF] for the latter. More recently, Yang et al. introduced a new semi-supervised constraint score which takes advantage of the local geometrical structure of unlabeled data samples as well as constraints deduced from prototypes [START_REF] Yang | Semi-supervised feature selection for audio classification based on constraint compensated laplacian score[END_REF][START_REF] Yang | Semi-supervised minimum redundancy maximum relevance feature selection for audio classification[END_REF].

Because the above-mentioned constraint scores are part of the filter approach, they all evaluate features one by one [START_REF] Sheikhpour | A survey on semi-supervised feature selection methods[END_REF]. The score of a feature subset is estimated by the sum of the individual feature scores, and the evaluation of a feature subspace ignores correlations between features. Thus, learning algorithms that operate in a subspace of individually relevant features do not necessarily provide favorable results [START_REF] Benabdeslem | Efficient semi-supervised feature selection: constraint, relevance, and redundancy[END_REF]. In addition, the constraint scores proposed in the literature are based on the Laplacian of a similarity matrix. Because the similarity matrix is computed in the original feature space, state-of-the-art feature scores can also be corrupted by the curse of dimensionality.

In this paper, we propose a new constraint score that evaluates the relevance of features in the context of both supervised and semi-supervised learning. Our score assesses the ability of features to respect the available set of pairwise constraints. As this score can be used for feature selection, it is based on a similarity matrix that is computed in the considered feature subspace. Unlike existing constraint scores that are applied to each feature, our score can evaluate a subset of several features simultaneously. The proposed score is then used as a criterion by a sequential forward selection scheme to identify the most relevant subset of features with tractable computation [START_REF] Siedlecki | On automatic feature selection[END_REF].

The performance of the constraint scores is measured by the classification accuracy of the test data commonly obtained by the nearest neighbor classifier.

Previous studies use the entire training dataset with true class labels as prototypes by the classifier, while only few prototypes are used by the constraint scores. By using conditions that are similar to those in real-life applications, in this paper we propose using only available information. In the supervised context, only the prototypes involved in pairwise constraint generation are used by the classifier. In the semi-supervised context, we follow the same strategy proposed by Kalakech et al. [START_REF] Kalakech | Constraint score evaluation for spectral feature selection[END_REF] that first uses the constrained K-means algorithm [START_REF] Wagstaff | Constrained k-means clustering with background knowledge[END_REF] to classify the unlabeled training data samples and then uses the classified samples as prototypes in classifying the test data. However, in our case, instead of using the constrained K-means algorithm, we use constrained spectral clustering, which is based on the same similarity matrix concept used by the constraint scores.

The remainder of this paper is organized as follows. Section 2 provides a brief definitions on spectral graph theory and pairwise constraints generation. In Section 3 a primary state-of-the-art constraint scores is presented. Our proposed constraint score and the feature selection procedure are presented in Section 4. Experimental results achieved with benchmark databases pertaining to supervised and semi-supervised feature selection are provided and discussed in Section 5.

Preliminaries

Constraint scores are based on the concepts of spectral graph theory and pairwise constraints. In this section, we briefly give some notations and definitions related to these two concepts.

Spectral graph theory

Let X = [x 1 , x 2 , . . . , x n ] ∈ R n×d denote the set of n training data samples defined in a d-dimensional feature space, where x i = [x i1 , x i2 , . . . , x ir , . . . , x id ] ∈ R d is the i-th data sample of X, and x ir ,(r = 1, . . . , d) is the r-th feature value of the i-th data sample. Let F d = {f 1 , f 2 , . . . , f r , . . . , f d } denote the set of d feature vectors of X, where f r = [x 1r , x 2r , . . . , x nr ] T ∈ R n is the r-th feature vector.

In spectral graph theory, dataset X is represented by a complete undirected weighted similarity graph G = (V, E, W ) in which each data sample in X corresponds to a node. V is a non-empty set that contains all nodes, and E is the set of edges between any two nodes in V . Each edge in E is weighted by a similarity value w ij (i, j = 1, 2, . . . , n) between two nodes. The similarity matrix W that gathers similarities between all pairs of nodes is positive semi-definite and symmetric. Generally, the similarity w ij between two data samples x i and x j is computed by the following Gaussian kernel function [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]:

w ij = exp - δ 2 (x i , x j ) 2σ 2 i, j = 1, 2, . . . , n (1)
where δ(x i , x j ) is the Euclidean distance between two data samples x i and x j , and σ is a scaling parameter.

Dataset X can also be represented by the similarity matrix W KNN of the nearest neighbor subgraph G KNN whose node v i is connected to v j when x j is one of the K-nearest neighbors (KN N s) of x i such that:

w KNN ij =    w ij if x i ∈ KNN(x j ) or x j ∈ KNN(x i ) 0 otherwise (2)
The constraint scores are often formulated from a given unnormalized Laplacian matrix of W , that is defined as follows:

L = D -W (3)
where D ∈ R n×n is the diagonal degree matrix whose elements are d ii = n j=1 w ij [START_REF] Luxburg | A tutorial on spectral clustering[END_REF][START_REF] Shi | Normalized cuts and image segmentation[END_REF].

Pairwise constraint

In supervised and semi-supervised learning frameworks, often a limited amount of information on the training dataset is available. This prior knowledge is expressed by only few labeled data samples. It can also be expressed by pairwise constraints that mention if data samples belong to the same class (must-link)

or to different classes (cannot-link). Pairwise constraints can be given by the user or are generated from labeled data samples.

In this paper, we consider that only few labeled data samples (prototypes) characterize the k classes ω l , l = 1, . . . , k. Let X l ∈ R p×d , (X l ⊂ X) be the set of p prototypes that are associated with class ω l . From the overall set of prototypes

denoted X P X P = l=1,...,k X l , we can build set M of (k•p 2 -k•p = k•p•(p-1))
must-link pairs that are composed of two prototypes belonging to the same class:

M = (x i , x j ) ∈ X 2 | ∃ l = 1, . . . , k so that x i ∈ X l and x j ∈ X l (4)
We can also build set C of (k • (k -1) • p 2 ) cannot-link pairs that are composed of two prototypes belonging to different classes:

C = (x i , x j ) ∈ X 2 | ∃ (l, m); l = m; so that x i ∈ X l and x j ∈ X m (5)
Of all possible data pairs that can be extracted from X, those belonging to M or C are called constrained pairs, while the remaining pairs are called unconstrained pairs. Furthermore, data that are not prototypes are called unlabeled data samples and are gathered in subset X U = X/X P .

In the context of spectral graph theory, two graphs G M and G C are built from the sets of must-link constraints M and cannot-link constraints C. The similarity matrices W M ∈ R n×n and W C ∈ R n×n are defined as follows:

w M ij =    1 if (x i , x j ) ∈ M 0 otherwise (6) w C ij =    1 if (x i , x j ) ∈ C 0 otherwise (7)

Constraint scores

The performance achieved by learning algorithms such as classification or clustering depends on similarity, which is based on the Euclidean distance in the original d-dimensional feature space. Because these features are not always relevant, many authors select the best ones based on constraint scores that combine the concepts of spectral graph theory and pairwise constraints.

Constraint scores for supervised feature selection

In a supervised learning framework, only pairwise constraints are used to evaluate the relevance of features while unconstrained data samples are ignored.

Two constraint scores C 1 r and C 2 r of feature f r are defined by Zhang et al. [START_REF] Zhang | Constraint score: A new filter method for feature selection with pairwise constraints[END_REF] as follows:

C 1 r = (xi,xj )∈M (x ir -x jr ) 2 (xi,xj )∈C (x ir -x jr ) 2 = n i=1 n j=1 (x ir -x jr ) 2 w M ij n i=1 n j=1 (x ir -x jr ) 2 w C ij (8) C 2 r = (xi,xj )∈M (x ir -x jr ) 2 -λ (xi,xj )∈C (x ir -x jr ) 2 = n i=1 n j=1 (x ir -x jr ) 2 w M ij -λ n i=1 n j=1 (x ir -x jr ) 2 w C ij ( 9
)
where λ ∈ [0, 1] is a regularization parameter. In our experiments, λ is set to 1, as in [START_REF] Zhang | Constraint score: A new filter method for feature selection with pairwise constraints[END_REF]. These scores assume that the distance between the must-link data samples should be as low as possible, while the distance between the cannot-link data samples should be as high as possible.

The two scores can be formulated from the unnormalized constrained Lapla-

cian matrices L M = D M -W M and L C = D C -W C : C 1 r = f T r L M f r f T r L C f r (10)
C 2 r = f T r L M f r -λf T r L C f r (11)
Lower scores indicate a more relevant feature.

Constraint scores for semi-supervised feature selection

Semi-supervised feature selection involves the analysis of both pairwise constraints and unlabeled data samples. It considers both the discriminating power of the pairwise constraints and the local properties of the unlabeled data samples.

Zhao et al. [START_REF] Zhao | Locality sensitive semi-supervised feature selection[END_REF] introduced the semi-supervised constraint score C 3 r , called the locality sensitive discriminant analysis score, which combines the similarity matrix W C constructed from the cannot-link constraints (Eq. ( 7)) and the similarity matrix W KNN1 ∈ R n×n which is constructed from the set of must-link constraints and unlabeled data samples as follows:

w KNN1 ij =          γ if (x i , x j ) ∈ M 1 if (x i ∈ X U or x j ∈ X U ) and x i ∈ KNN(x j ) or x j ∈ KNN(x i ) 0 otherwise ( 12
)
where γ is a constant parameter. In our experiments, γ is set to 100 and K is set to 5 as in [START_REF] Zhao | Locality sensitive semi-supervised feature selection[END_REF].

The constraint score C 3 r is defined as

C 3 r = n i=1 n j=1 (x ir -x jr ) 2 w KNN1 ij n i=1 n j=1 (x ir -x jr ) 2 w C ij = f T r L KNN1 f r f T r L C f r (13)
where L KNN1 = D KNN1 -W KNN1 is the unnormalized Laplacian matrix of W KNN1 , D KNN1 being the degree matrix computed from W KNN1 .

This score ensures that two data samples related by a cannot-link constraint are well separated. It also implicitly takes into account the unlabeled data samples but favors pairs of must-link data samples by assigning them high weights in the matrix W KNN1 . Moreover, the similarity matrix W KNN1 represents the links between the KN N s of the unlabeled data samples by binary weighting them.

Kalakech et al. [START_REF] Kalakech | Constraint scores for semisupervised feature selection: A comparative study[END_REF] proposed a semi-supervised constraint score C 4 r that is less sensitive to the constraint sets by a simple combination of scores computed on prototypes and unlabeled data samples. More precisely, C 4 r attempts to identify a trade-off between the unsupervised Laplacian score L r and the supervised constraint score C 1 r (see Eq. ( 8)) by multiplying both scores :

C 4 r = L r • C 1 r ( 14
)
The Laplacian score L r of a feature f r is defined as in [START_REF] He | Laplacian score for feature selection[END_REF] as follows: L r favors features with high variance and tends to select those with a strong ability to preserve locality, while C 1 r seeks to select features with a strong ability to preserve pairwise constraints.

L r = n i=1 n j=1 (x ir -x jr ) 2 w ij n i=1 (x ir -fr ) 2 d ii = f T r L fr f T r D fr
Finally, C 4 r is defined in terms of Laplacian matrices as follows:

C 4 r = f T r L fr f T r D fr • f T r L M f r f T r L C f r (16)
L and D are deduced from similarity matrix W (Eq. ( 1)).

Benabdeslem and Hindawi proposed another constraint score called the constrained Laplacian score that combines the similarity matrices W C and W KNN2

for semi-supervised feature selection [START_REF] Benabdeslem | Constrained laplacian score for semisupervised feature selection[END_REF]. The similarity matrix W KNN2 is expressed as

w KNN2 ij =    w ij if (x i , x j ) ∈ M or x i ∈ KNN(x j ) or x j ∈ KNN(x i ) 0 otherwise (17)
The score C 5 r is defined as follows [START_REF] Benabdeslem | Constrained laplacian score for semisupervised feature selection[END_REF][START_REF] Benabdeslem | Ensemble constrained laplacian score for efficient and robust semi-supervised feature selection[END_REF]:

C 5 r = n i=1 n j=1 (x ir -x jr ) 2 w KNN2 ij n i=1 n j=1|∃l,(x l ,xj )∈C (x ir -α i jr ) 2 d KNN2 ii ( 18
)
where

α i jr =          x ir if (x i , x j ) ∈ C fr if i = j and x i ∈ X U x ir otherwise (19)
C 5 r represents an enhanced version of both scores: the Laplacian score L r and the supervised constraint score C 1 r . In fact, the Laplacian score L r can be seen as a special version of C 5 r when there are no labeled data samples (α i jr = fr ), and when (α i jr = x ir ), C 5 r can be considered as an adjusted version of the constraint score C 1 r . Subsequently, C 5 r is defined in terms of Laplacian matrices as follows [START_REF] Benabdeslem | Constrained laplacian score for semisupervised feature selection[END_REF]:

C 5 r = f T r L KNN2 f r f T r L C D KNN2 f r (20)
Finally, on the one hand, this score ensures that data samples that are neighbors or related by must-link constraints should be close together when they are projected to relevant features. On the other hand, the distance between data samples that are related by cannot-link constraints should be as high as possible.

A second semi-supervised constrained Laplacian score, referred to C 6 r , has been proposed by Benabdeslem and Hindawi in [START_REF] Benabdeslem | Efficient semi-supervised feature selection: constraint, relevance, and redundancy[END_REF] as follows:

C 6 r = n i=1 n j=1 (x ir -x jr ) 2 (w KNN ij + N ij ) n i=1 n j=1 (x ir -α i jr ) 2 d KNN ii = f T r L KNN3 f r f T r L C D KNN f r (21)
where the diagonal matrix D KNN is deduced from the similarity matrix W KNN (see Eq. 2). N ij is given as follows:

N ij =                -w ij if (x i , x j ) ∈ M and x i ∈ KNN(x j ) or x j ∈ KNN(x i ) w 2 ij if (x i , x j ) ∈ M and x i / ∈ KNN(x j ) and x j / ∈ KNN(x i ) or (x i , x j ) ∈ C and x i ∈ KNN(x j ) or x j ∈ KNN(x i ) 0 otherwise (22)
The Laplacian matrix L KNN3 is computed from the similarity matrix W KNN3 = w KNN ij + N ij defined as follows:

w KNN3 ij =                w 2 ij + w ij if (x i , x j ) ∈ C and x i ∈ KNN(x j ) or x j ∈ KNN(x i ) w 2 ij if (x i , x j ) ∈ M and x i / ∈ KNN(x j ) and x j / ∈ KNN(x i ) w ij if x i ∈ X U or x j ∈ X U and x i ∈ KNN(x j ) or x j ∈ KNN(x i ) 0 otherwise (23)
C 6 r combines the power of the local geometric structure offered by unlabeled data samples, with the constraint preserving ability offered by prototypes. The additional weight w 2 ij was introduced in order to more differentiate the features in the both bad cases, i.e., when two data samples are related by a must-link constraint but are not neighbors and when two neighboring data samples are related by a cannot-link constraint.

C 7
r , called constraint compensated Laplacian score, which takes advantage of the local geometrical structure of unlabeled data samples as well as constraint information deduced from labeled data samples [START_REF] Yang | Semi-supervised feature selection for audio classification based on constraint compensated laplacian score[END_REF][START_REF] Yang | Semi-supervised minimum redundancy maximum relevance feature selection for audio classification[END_REF]:

C 7 r = n i=1 n j=1 (x ir -x jr ) 2 (w KNN ij + Nij ) r + b r - w r (24)
Nij is given as follows:

Nij =                1 -w ij if (x i , x j ) ∈ M and x i ∈ KNN(x j ) or x j ∈ KNN(x i ) γw ij if (x i , x j ) ∈ C and x i ∈ KNN(x j ) or x j ∈ KNN(x i ) λ if (x i , x j ) ∈ M and x i / ∈ KNN(x j ) and x j / ∈ KNN(x i ) 0 otherwise (25)
where γ and λ are the parameters set to the empirical values of 0.9 and 0.5, respectively [START_REF] Yang | Semi-supervised feature selection for audio classification based on constraint compensated laplacian score[END_REF][START_REF] Yang | Semi-supervised minimum redundancy maximum relevance feature selection for audio classification[END_REF].

= n i=1 (x ir -fr ) 2 d ii (26) b r = k l=1 |X l |( f (l) r -f P r ) 2 (27) w r = k l=1 |X l |(σ (l) r ) 2 (28)
where

|X l | = p is the number of prototypes of the lth class, f P r = n i=1|xi∈X P xir k.p
is the mean of the rth feature of the labeled dataset,

f (l) r = n i=1|xi∈X l xir p and σ (l)
r denote the mean and variance of the rth feature of the lth class, respectively.

The semi-supervised constraint score C 7 r can be expressed in terms of Laplacian matrices as follows [START_REF] Yang | Semi-supervised minimum redundancy maximum relevance feature selection for audio classification[END_REF]:

C 7 r = 2(f r ) T L KNN4 f r ( fr ) T D KNN4 fr + 2( f P r ) T L P f P r -( f P r ) T D P f P r (29)
The Laplacian matrix L KNN4 is computed from the similarity matrix W KNN4 = W KNN + Nij that is expressed as

w KNN4 ij =                      1 if (x i , x j ) ∈ M and x i ∈ KNN(x j ) or x j ∈ KNN(x i ) λ if (x i , x j ) ∈ M and x i / ∈ KNN(x j ) and x j / ∈ KNN(x i ) (1 -γ)w ij if (x i , x j ) ∈ C and x i ∈ KNN(x j ) or x j ∈ KNN(x i ) w ij if x i ∈ X U or x j ∈ X U and x i ∈ KNN(x j ) or x j ∈ KNN(x i ) 0 otherwise (30)
L P = D P -W P is the unnormalized Laplacian matrix of W P , D P being the degree matrix computed from W P . The cells of W P are set to 1/|X l | when two data samples are prototypes that belong to the same class and to 0, otherwise:

w P ij =    1/|X l | if x i ∈ X l and x j ∈ X l 0 otherwise (31) Note that f P r = f P r -f P r .
The features providing the lowest scores C 3 , C 4 , C 5 C 6 and C 7 are the most relevant features. However, because these scores are based on the Laplacian matrices and the diagonal degree matrices, which are deduced from the similarity matrices, they are computed in the original d-dimensional feature space and can be corrupted by the curse of dimensionality.

Proposed constrained feature selection

Existing constraint scores estimate the relevance of each feature considered independently and separately from each other. Because these scores do not take into account the correlation between features, we propose a new constraint score that estimates the relevance of a subset of features at once.

Proposed constraint score

Unlike previous constraint scores that are based on Laplacian matrices, we propose a constraint score that is based only on similarity matrices to evaluate the relevance of a subset of m features denoted

F m = {f 1 , . . . , f m } (m = 1, 2, . . . , d).
The proposed score, denoted ε * (F m ), can be used in one of two learning contexts: * = supervised (S) or semi-supervised (SS). In the supervised learning context, we use only the must-link and cannot-link constraint sets to select features using ε S (F m ). In the semi-supervised learning context, we exploit both the pairwise constraint sets and the information contributed by the unlabeled data samples to select features using ε SS (F m ).

The relevance of F m is evaluated by means of the distance between the target similarity matrix Ŵ * , * = S or SS, which is defined from the given constraints and a similarity matrix W (F m ) that is computed with the subset F m of features. The score ε * (F m ), which should be as low as possible, is expressed as the following square error:

ε * (F m ) = W (F m ) -Ŵ * 2 with * = S or SS (32)
where • 2 is the Euclidean norm. Thus, ε * (F m ) can be rewritten as

ε * (F m ) = n i=1 n j=1 w ij (F m ) -ŵ * ij 2 (33)
W (F m ) ∈ R n×n is the similarity matrix computed on the dataset X with the subset of features F m :

w ij (F m ) = exp - δ 2 x (m) i , x (m) j 2σ 2 i, j = 1, 2, . . . , n (34) where x
(m) i is the vector of the i-th data point characterized by the subset F m .

Ŵ * ∈ R n×n is the target matrix whose cells correspond to must-link pairs are set to 1, while the cells corresponding to the cannot-link pairs are set to 0.

Supervised constraint score

For supervised learning, the target similarity matrix Ŵ S is defined as follows:

ŵS ij =          1 if (x i , x j ) ∈ M 0 if (x i , x j ) ∈ C w ij (F m ) otherwise (35)
The cells of Ŵ S corresponding to the unconstrained pairs are set to w ij (F m ) so that they are not taken into account by ε S (F m ), (w ij (F m )-ŵS ij = 0). Therefore, from Eqs. (33) and ( 35), it can be easily demonstrated that ε S (F m ) can be expressed as

ε S (F m ) = (x i , x j ) ∈ M (x i , x j ) ∈ C w ij (F m ) -ŵS ij 2 (36)
Finally, ε S (F m ) makes it possible to select the set of features with the best constraint preserving ability and without taking into account the unlabeled data samples. The concept behind ε S (F m ) is simple and natural. A relevant feature subset produces similarity w ij (F m ) between two must-link samples that is close to 1; however, a good feature subset should provide a similarity between cannot-link samples close to 0.

Semi-supervised constraint score

For semi-supervised learning new must-link pairs are constructed from the prototypes and unlabeled data samples to compute the binary target similarity matrix Ŵ SS :

ŵSS ij =    1 if (x i , x j ) ∈ M SS 0 otherwise (37)
The cell (i, j) of Ŵ SS is set to 1 (0 otherwise) when (x i , x j ) belongs to M SS , which is a new set of must-link pairs deduced from prototype subsets X l , l = 1, . . . , k and unlabeled data samples as follows:

M SS = (x i , x j ) ∈ X 2 |∃l = 1, . . . , k so that N P (x i ) ∈ X l and N P (x j ) ∈ X l (38)
The nearest prototype N P (x i ) is the prototype whose distance from a sample data point x i ∈ X in the original d-dimensional feature space is the smallest one:

N P (x i ) = arg min y ∈ l=1,...,k X l δ 2 (x i , y) (39)
Thus, pair (x i , x j ) belongs to M SS when the nearest prototype of x i and x j both belong to the same subset X l of prototypes. Because N P (x i ) is x i when

x i belongs to prototype subset X l , set M is included in M SS (M ⊂ M SS ).
Therefore, M SS increases the contribution of the pairwise constraints M for semi-supervised feature selection. An illustration of M SS is provided in Fig. 1, of Eq. ( 12), W KNN2 of Eq. ( 17), W KNN3 of Eq. ( 23) and W KNN4 of Eq. ( 30), which is used by the semi-supervised constraint score C 3 , C 5 , C 6 and C 7 , respectively. These links are established from the pairwise constraints and the K-nearest neighbors (KN N s) of the data samples. To compare them with our set M SS , we set K to 1 for W KNN1 , W KNN2 , W KNN3 and W KNN4 and

in
define 1NN = {(x i , x j ) ∈ X 2 | x i ∈ X U or x j ∈ X U but x i = 1NN(x j ) or x j =
1NN(x i )} to describe links that are generated by these four matrices. Figure 1 demonstrates that M SS better respects the geometric structure of classes than evaluate the relevance of each feature, the selection procedure ranks features with respect to one of the scores to identify the most relevant sets of features.

W KNN1 , W KNN2 , W KNN3
For illustration purposes, we consider the Dermatology database from UCI repository [START_REF] Lichman | UCI machine learning repository[END_REF], which contains 366 samples characterized by 34 features and is regrouped into 6 classes. We randomly pick out p = 3 prototypes for each class. Because our score ε * can evaluate the relevance of a feature subset at once, it is more advantageous to use a selection procedure capable of combining features between them. For this, we use the sequential forward feature selection technique due to its simplicity [START_REF] Siedlecki | On automatic feature selection[END_REF].

To evaluate the relative relevance of d features, we first consider each feature uates the features jointly and considers the dependency among them by virtue of our proposed constraint score, it belongs to the filter approach and retains all the advantages of this approach comparatively to wrapper approach (i.e., independence from the classifier, high computation speed, simple and rapid implementation).

Experiments on benchmark databases

We evaluated and compared our proposed constraint score with several constraint scores and several well-known feature selection methods on datasets originating from benchmark databases. We first examined the supervised feature scores (ε S , C 1 , and C 2 ). Then, we assessed the performance attained by semisupervised feature scores (ε SS , C 3 , C 4 , C 5 , C 6 and C 7 ). Because the data were scaled between 0 and 1, the scaling parameter σ used to compute the similarity matrices was set to 1. Feature selection procedures were performed on the training datasets and repeated over 100 runs. In each feature selection run, we automatically generated a set of pairwise constraints as follows. For each class, we randomly selected from the training dataset X, k prototype subsets X l (|X l | = p). Then, we deduced sets M , C, and M SS of the pairwise constraints using Eqs. ( 4), [START_REF] Zhang | Constraint score: A new filter method for feature selection with pairwise constraints[END_REF], and (38), respectively.

Databases

Our experiments were performed using twelve six well-known and commonly used benchmark databases in the feature selection framework, namely, the Wisconsin Breast Cancer (WBCD), Image Segmentation, Wisconsin Diagnostic Breast Cancer (WDBC), Ionosphere, Dermatology, Libras Movement, Multiple

Features and CNAE-9 databases from the UCI repository [START_REF] Lichman | UCI machine learning repository[END_REF], and the Olivetti Research Laboratory (ORL), Yale, Pie10P, and ALLAML databases from the ASU feature selection repository [START_REF] Zhao | Advancing feature selection research[END_REF]. The details of these databases are provided in Table 1. All the databases had numerical features, and the class label of each data sample was clearly defined. We normalized each dataset between 0 and 1 to ensure that the scales of all features were equal. For WBCD, 410 samples were used as the training dataset and 273 samples were used as the test dataset, while for WDBC, 376 samples were used as the training dataset and 193 samples were used as the test dataset [START_REF] Sheikhpour | A survey on semi-supervised feature selection methods[END_REF]. For the remaining databases, Ionosphere, Dermatology, Libras Movement, and ORL, we used half of the data samples from each class as the training dataset and the remaining data samples as the test dataset [START_REF] Kalakech | Constraint scores for semisupervised feature selection: A comparative study[END_REF]. It should be noted that the classes were not equiprobable for all databases with the exception of the WBCD, WDBC, Ionosphere, Dermatology and ALLAML Libras Movement and ORL databases. We extracted (k • p) prototypes from the training sample set of each database with p ranging from 2 to 4. Table 1 also displays the number (k • p) of prototypes extracted from the training sample set of each database with p ranging from 2 to 4. These prototypes were used to build pairwise constraints and select relevant features.

They were also used in the classification step to evaluate the relevance of the selected features. It can be seen that the curse of dimensionality is faced, as the number of prototypes always remains less than the number of features d.

Supervised feature selection

For supervised learning, only k prototype subsets X l are used to select features. The performance obtained by the proposed supervised constraint score ε S is compared first with that attained by supervised constraint scores C 1 and ReliefF [START_REF] Robnik-Šikonja | Theoretical and empirical analysis of relieff and rrelieff[END_REF] and minimal-redundancy-maximal relevance (mRMR) [START_REF] Peng | Feature selection based on mutual information criteria of max-dependency, max-relevance, and min-redundancy[END_REF]. In addition to these filter feature selection methods, a wrapper-type supervised feature selection [START_REF] Li | Using the k-nearest neighbor algorithm for the classification of lymph node metastasis in gastric cancer[END_REF] is used for the comparison.

The relevance of the selected subset of features is evaluated by the accuracy measure of the test dataset. Each test data point is projected onto the retained feature space and is assigned to one of the k classes according to the nearest neighbor rule. In general, the entire training dataset is used as prototypes by the nearest neighbor classifier to classify the test data, whereas only few prototypes are used by the supervised constraint scores [START_REF] Kalakech | Constraint scores for semisupervised feature selection: A comparative study[END_REF][START_REF] Zhang | Constraint score: A new filter method for feature selection with pairwise constraints[END_REF][START_REF] Benabdeslem | Constrained laplacian score for semisupervised feature selection[END_REF]. Here, we propose performing the selection and evaluation with only the same available k prototype subsets X l (|X l | = p). In this way, these two steps are performed in conditions similar to those of real-life applications.

Because the accuracy highly depends on the number of selected features and prototypes, we first evaluate accuracy with respect to the number of features m for a given number of prototypes p. Then, for a relevant number of features, we compare the accuracy obtained by the supervised constraint scores with different numbers of prototypes.

Accuracy versus number of selected features

Figure 4 illustrates the average accuracy for ε S , C 1 , and C 2 obtained when the number of prototypes p is set to 3. From this figure, it can be seen that in most cases, the average accuracy obtained with the proposed scores ε S is higher than that obtained with constraint scores C 1 and C 2 . The curves of ε S differ from the curves of C 1 and C 2 for the Image Segmentation, Ionosphere, Dermatology, and Libras Movement databases. However, the curves of ε S , C 1 , and C 2 overlap for the other databases. Because the performance of these scores is averaged over 100 runs with different generations of constraints, comparing them is difficult. As a result, we compare these scores by examining their accuracy at each of the 100 runs.

For a fixed number of selected features, in each of the 100 runs, we propose ranking the supervised constraint scores ε S , C 1 , and C 2 in descending order of accuracy. Let us denote rank

[c]
q the rank of the supervised constraint scores c = ε S , C 1 , or C 2 at run q. The score with the highest accuracy is ranked 1, while the score with the lowest accuracy is ranked 3. Scores with the same accuracy have the same rank. We compute the rank sum R [c] for each score as follows:

R [c] = 100 q=1 rank [c] q (40)
The method with the lowest rank sum is considered to be the score that provides the best results.

For the WBCD, Image Segmentation, WDBC, Ionosphere, Dermatology, Libras Movement, Multiple Features, and ORL databases, the accuracy of the supervised constraint scores appears to be stable when the number of selected features is higher than 5, 11, 12, 13, 15, 45, 100, and 300, respectively (see Fig. 4).

The rank sum of each supervised constraint score is then computed by considering the first 5, 11, 12, 13, 15, 45, 100, and 300 features for the WBCD, Image In Table 2, the row with an asterisk for each of the eight six databases, provides the rank sum when the number of prototypes p is set to 3. These rows indicate that ε S provides the lowest rank sum for all databases (indicated in bold) except for the WDBC database. Thus, the accuracy provided by the features selected by ε S is higher than those obtained by C 1 and C 2 .

Accuracy versus number of prototypes

In this section, we compare the accuracy obtained by the supervised constraint scores ε S , C 1 , and C 2 for various number of prototypes p. For this purpose, we use the number of selected features that is provided in the previous section. Because the number of prototypes p must be higher than or equal to 2

to generate at least one must-link constraint by class, it ranges from 2 to 4.

Table 2 displays the rank sums obtained for the eight databases and demonstrates that the features selected by ε S provide higher accuracy rates than those obtained by the features selected by C 1 and C 2 . In fact, the score ε S provides the lowest rank sum (indicated in bold) 19 times out of the 24 rows in Table 2.

The improvement provided by ε S has two main causes. First, because ε S estimates the relevance of a subset of features, it takes into account the correlation between them, whereas C 1 and C 2 ignore this correlation. Second, ε S computes the similarity matrix between samples in the considered subset of features (see.

Eq. ( 36)) whereas C 1 and C 2 do not compute the similarity between samples.

Comparison with other supervised feature selection methods

To further illustrate the effectiveness of our supervised constraint score ε S , it is compared with well-established supervised feature selection methods including ReliefF [START_REF] Robnik-Šikonja | Theoretical and empirical analysis of relieff and rrelieff[END_REF] and mRMR [START_REF] Peng | Feature selection based on mutual information criteria of max-dependency, max-relevance, and min-redundancy[END_REF]. A wrapper-type supervised feature selection that selects features in a forward sequential manner by means of nearest neighbor classifier (SFS-1NN) is also used for the comparison [START_REF] Li | Using the k-nearest neighbor algorithm for the classification of lymph node metastasis in gastric cancer[END_REF]. Figure 5 displays the average accuracies over 100 runs achieved by the supervised feature selection methods on high dimensionality databases compared with the r% selected features, when the number p of prototypes is set to 3. r is set to 5 for datasets with a number d of features greater than 3000, 10 for datasets with a number of features in the range [2000 -3000], 20 for datasets with a dimensionality in the range [300 -2000] and 100 for lower dimensionality (less than 300). Figure 5 shows that the averaged accuracy on the features selected by ε S is usually higher than that of the filter methods but remains lower than that of the wrapper method SFS-1NN. This observation is confirmed by the averaged accuracies over the top r% features that are reported in Table 3. To better assess the results obtained for each algorithm, we ranks the algorithms for each dataset separately, the best performing algorithm obtaining the rank of 1, the second best rank 2, etc. In case of ties, we assigns the same rank. Finally, the average rank over all databases is depicted in the bottom row of Table 3. We can find overall that our score ε S outperforms the four other filter methods, but remains behind the wrapper method.

Semi-supervised feature selection

Semi-supervised learning analyzes the k prototype subsets X l in addition to unlabeled samples that belong to the learning dataset X. Because our score SS uses unlabeled samples to construct the new set M SS of must-link pairwise constraints, it is interesting to first evaluate the relevance of M SS . Then, the performance obtained by the proposed semi-supervised constraint score ε SS is compared with that attained by the semi-supervised constraint scores C c (c = 3, . . . , 7) and by two semi-supervised feature selection methods namely the semisupervised pairwise constraint-guided sparse (SCGS) learning method [START_REF] Liu | Pairwise constraint-guided sparse learning for feature selection[END_REF] and the ensemble constrained Laplacian score (EnsCLS) method [START_REF] Benabdeslem | Ensemble constrained laplacian score for efficient and robust semi-supervised feature selection[END_REF].

For this purpose, we follow the same experimental scheme as that of supervised learning. First, we evaluate the accuracy with respect to the number of features m for a given number of prototypes p. Then, for a number of features m that are considered relevant, we compare the accuracy achieved by the semi-supervised constraint scores with different numbers of prototypes.

Correct must-link pairs of M SS

The proposed semi-supervised constraint score ε SS builds the set of new must-link pairs M SS (see Eq. ( 38)), which extends set M . To assess the resulting set M SS , we propose comparing it with M by computing the rate of new mustlink pairs (N M L). Then, the rate of correct new must-link pairs (CN M L) is computed as follows to evaluate the relevance of these new constraints:

N M L = |M SS | |M | (41) CN M L = | {(x i , x j ) ∈ M SS \ M | ω(x i ) = ω(x j )} | |M SS \ M | (42)
where ω(x i ) is the true class of the data sample x i . Here, N M L ≥ 1 and 0 ≤ CN M L ≤ 1, and higher N M L and CN M L lead to an improved M SS . Therefore, we can deduce that the majority of these new must-link pairwise constraints should be used to select features.

Accuracy versus number of selected features

The performance of semi-supervised constraint scores is measured according to the classification accuracy of the test dataset. In semi-supervised learning, the evaluation step is generally performed in the supervised context, in which the training dataset is used as prototypes by the nearest neighbor classifier to classify the test data [START_REF] Kalakech | Constraint scores for semisupervised feature selection: A comparative study[END_REF][START_REF] Zhao | Locality sensitive semi-supervised feature selection[END_REF][START_REF] Benabdeslem | Constrained laplacian score for semisupervised feature selection[END_REF]. Instead, the features are selected in the semisupervised learning context in which only few prototypes are considered. To perform selection and evaluation in the same semi-supervised context that is similar to real-life applications, we follow the strategy proposed by Kalakech et al. [START_REF] Kalakech | Constraint score evaluation for spectral feature selection[END_REF] by applying the constrained spectral clustering detailed in [START_REF] Kamvar | Spectral learning[END_REF] instead of the constrained K-means algorithm to classify the unlabeled training data samples. Then, we use the sample as prototypes to classify the test data subset.

Finally, each test data point is projected onto the same retained feature subspace and is assigned to a class using the nearest neighbor rule, which uses the available prototypes in addition to the training data samples that have been previously classified.

Figure 6 presents the average accuracy for ε SS , C 3 , C 4 , C 5 , C 6 and C 7 obtained in the eight test sets using the semi-supervised learning evaluation over 100 runs when the number of prototypes p was set to 3. This figure indicates that our semi-supervised constraint score ε SS provides higher accuracy than that achieved by the other semi-supervised constraint scores C c (c = 3, . . . , 7).

The curves of ε SS clearly differ from those of C c (c = 3, . . . , 7), which overlap.

We also determined the rank sum R [c] obtained by the semi-supervised constraint scores ε SS and C c (c = 3, . . . , 7) for the eight databases over 100 runs and ORL databases, respectively (see Fig. 6). The rows with an asterisk in Table 5 indicate that ε SS provides the lowest rank sum for all the databases (indicated in bold) except for the WDBC Ionosphere database. Thus, the accuracy provided by the features selected by ε SS is higher than that obtained by C c (c = 3, . . . , 7) in the majority of cases.

Accuracy versus number of prototypes

In this section, we compare the accuracy obtained by the semi-supervised constraint scores ε SS , C 3 , C 4 , C 5 , C 6 and C 7 versus the number of prototypes p.

For this purpose, we use the number of selected features provided in Sect.

Table 5 displays the rank sum over 100 different sets of prototypes when the number of prototypes p ranges from 2 to 4. It can be seen that ε SS provides the lowest rank sum 20 out of 24 times. These results indicate that the accuracy provided by features selected by ε SS is higher than the accuracy of those obtained by C c (c = 3, . . . , 7) in most cases.

In addition, ε SS outperforms C c (c = 3, . . . , 7) because it takes into account the correlation between features. Furthermore, it computes the similarity matrix between samples in the selected feature subspace, while C c (c = 3, . . . , 7) compute the similarities in the high-dimensional original feature space.

Comparison with other semi-supervised feature selection methods

To further illustrate the effectiveness of our semi-supervised constraint score ε SS , it is compared with other semi-supervised feature selection methods, namely the SCGS [START_REF] Liu | Pairwise constraint-guided sparse learning for feature selection[END_REF] and EnsCLS [START_REF] Benabdeslem | Ensemble constrained laplacian score for efficient and robust semi-supervised feature selection[END_REF]. In addition, two unsupervised feature selection methods including the Laplacian score (LS) [START_REF] He | Laplacian score for feature selection[END_REF] and the spectral feature selection (Spec) method [START_REF] Zhao | Spectral feature selection for supervised and unsupervised learning[END_REF], which are based on graph theory and similarity matrices, are used for the comparison. Note that EnsCLS combines both a over the top r% features for all databases, are reported in Table 6. The averaged rank of each algorithm over all datasets, is reported in the bottom row of the Table 6. From Fig. 7 and Table 6, we can make the following observations:

-When we compare ε SS with semi-supervised features selection methods, we notice that SCGS outperforms ε SS for only 3 out 12 databases (CNAE9, Pie10P

and ALLAML), and EnsCLS works better than ε SS for only two databases (Ionosphere and Pie10P). ε SS outperforms the two constraint scores C 4 and C 7 for all databases except for the ALLAML database where C 7 appears more efficient.

-By comparing ε SS with the unsupervised features selection methods, we find that ε SS outperforms LS and Spec for all databases except for ALLAML where Spec produces a better averaged accuracy and WDBC where LS displays the same averaged accuracy than ε SS .

However, overall, ε SS outperforms all of these competitive methods given its lower average rank value. In some cases, ε SS performs better results than the wrapper method SFS-1NN as for Multiple Features, Yale and ORL databases.

Discussion

In this section, we first discuss the influence of the scaling parameter σ, the time running and the optimal number of selected features of the proposed constraint scores ε S and ε SS for feature selection. Then, we compare the accuracy provided by the supervised constraint score ε S and the semi-supervised constraint score ε SS .

Influence of the scaling parameter σ

Like the constraint scores C c , (c = 1, 2, ..., 7), our constraint scores ε S and ε SS depend on the scaling parameter σ used to compute the similarity matrices (see Eqs. 32 to 34). In previous experiments, the scaling parameter σ was set to 1. To analyze its effect on the performance of our constraint scores, we display in Fig. 8 the average accuracy over 100 runs obtained by our constraint scores ε S and ε SS on the both WBCD, WDBC and Dermatology databases with a number of prototypes p set to 3, when σ ranges from 0.2 to 1.4. To avoid the influence of the spectral clustering, the accuracy is assessed on the test dataset under the same conditions for the both scores, i.e., by using the same nearest neighbor classifier based only on the k prototype subsets X l (|X l | = p). The curves show that overall, σ has a significant influence on the results, but that a value between 0.6 and 1.2 is best suited to the experimented databases.

Time running

Our proposed score ε * has the advantage to evaluate the relevance of a feature subset at once thanks to the sequential forward feature selection technique.

In order to quantify the computational effort required by the feature selection based on our constraint scores ε S and ε SS , its running time is compared with that consumed by the wrapper method SFS-1NN that also uses the sequential forward selection as search algorithm. The running times of the mRMR and EnsCLS methods are also represented as an indication. Table 7 displays the computation times obtained by the above algorithms on some datasets of Table 1 with a Intel Core i7 3.60GHz, with 8GB RAM computer. Although the running times are computer dependent, they give an idea of the computation time required by the algorithms for various dimensionalities d, class numbers k, sample numbers n and numbers m of selected features. ε S and ε SS are relatively more time-consuming than mRMR and EnsCLS, even if that of ε S remains comparable. In fact, the computation time of ε S and ε SS in itself is not very high, but it is rather the sequential forward selection procedure which increases their computation time. For ε SS , the high computational time is mainly due to the procedure for propagating labels on unlabeled data samples, the complexity of which is proportional to the size of the data samples. Indeed, we can see on Table 7 that the time consumed by ε SS on the Multiple Features dataset which contains 2000 data samples of dimension 649 is much larger than that obtained on the ALLAML dataset which contains 72 data samples of higher dimension (7129). However, ε S and ε SS remain always faster than the wrapper method SFS-1NN. This proves that the calculation of scores is faster than the calculation of the accuracy by the 1NN classifier. It should be emphasized that the computational time cannot be considered as a crucial drawback because it is always possible to reduce it by adopting a faster search strategy, instead of the sequential forward selection, or by reducing the number of unlabeled data samples by clustering as in [START_REF] Benabdeslem | Efficient semi-supervised feature selection: constraint, relevance, and redundancy[END_REF], [START_REF] Benabdeslem | Ensemble constrained laplacian score for efficient and robust semi-supervised feature selection[END_REF]. to automatically determine the optimal number of features. Indeed, as mentioned in Sec. 4.2, the curve of constraint scores ε S and ε SS versus number of features presents a minimum that can be considered as the optimal number of features. To illustrate this strategy, in Fig. 9 we simultaneously display the average curves of ε S and ε SS and the average accuracies obtained with ε S and ε SS versus number of features on Dermatology database. We can show that the average accuracy Acc( m) obtained with the subset of features F ( m) that correspond to the minimum of ε SS coincides or is close to the maximum value of accuracy Acc(m * ). The average accuracies Acc( m) and Acc(m * ) as well as the corresponding optimal number m and m * of features obtained on all databases are depicted in Table 8. From Table 8, we find that, in most cases, the average accuracies Acc( m) and Acc(m * ) are close but with a number m of features less than m * . 4 indicates that CN M L is close to 0.7 for all databases except for the Image Segmentation, Ionosphere and Libras Movement databases, in which CN M L remains less than or equal to 0.6. It can be deduced from this table that the improvement caused by ε SS compared with ε S is dependent on the relevance of M SS compared to that of M . When most of the new must-link constraints in M SS are correct, our semi-supervised constraint score outperforms the supervised constraint score in the majority of cases.

Conclusion

In this paper, we presented a new constraint score for feature selection in the context of both supervised and semi-supervised learning. This score evaluates a subset of features at one time, whereas state-of-the-art constraint scores evaluate only one feature at a time. This makes it possible to identify redun- A possible way to increase the performance of our scores in terms of accuracy is to adaptively determine the scaling parameter σ. One solution envisaged is to proceed as for the calculation of our score, i.e., by means of the distance between the target similarity matrix defined from the given constraints and a similarity matrix that is computed with the candidate values of σ. Another possible solution for improving the relevance of selected features is to ensure of the coherence of pairwise constraints.

The proposed scores, especially in the semi-supervised context, are timeconsuming. However, the computational time cannot be considered as a crucial drawback because it is always possible to reduce it. We plan to further investigate this in our future work.

In this paper, we consider that prior knowledge is represented by class prototypes from which pairwise constraints are deduced. In this case, we can compute new pairwise constraints from prototypes, and all available constraints contribute to efficient feature selection. In future work, we intend to generalize our score to prior knowledge that is formalized only by pairwise constraints. This

1 T

 1 D1 , fr = f r -fr and 1 = [1, . . . , 1] T .

r

  is the variance of the rth feature of the whole dataset X, variance and between-class variance of the rth feature of the labeled dataset (protoypes) X P , respectively.

  r

  which only must-link pairs are represented for clarity. For comparison purposes, Fig. 1(d) displays the links corresponding to similarity matrix W KNN1

  and W KNN4 1NN by taking into account the margin between unlabeled samples and prototypes of different classes.

Figure 1 :

 1 Figure 1: Must-link pairs in a semi-supervised learning context. (a) Set of unlabeled data samples (displayed with color of class labels for a better interpretation of results) and two prototype sets X 1 and X 2 . (b) Set M . (c) Set M SS . (d) Set 1NN Links corresponding to similarity matrices W KNN1 , W KNN2 , W KNN3 and W KNN4 .

Finally, ε

 ε SS (F m ) = n i=1 n j=1 w ij (F m ) -ŵSS ij 2makes it possible to assess the ability of a feature subset to preserve the pairwise constraints provided by the user and the extended pairwise constraints provided by M SS . These extended pairwise constraints well represent the geometric structure of the classes.The supervised constraint scores C 1 and C 2 do not compute the similarity between prototypes in any feature space, while the classical semi-supervised constraint scores C c (c = 3, . . . , 7) evaluate the similarity between samples only in the original feature space. In contrast, our scores ε S and ε SS evaluate the similarity matrix W (F m ) in the selected feature subspace F m . This difference is essential because the classification or clustering of the data samples is performed in the selected feature subspace.4.2. Feature selection procedureFeature selection reduces data dimensionality by selecting the most relevant sets of original features. Because the state-of-the-art scores C c (c = 1, . . . , 7)

Figure 2

 2 Figure 2 presents the variation in scores C c (c = 1, . . . , 7) with respect to the number m of features. It should be noted that the score of a subset composed of the best m features in Fig. 2 is computed as the cumulative sum of the lowest scores divided by m. It can be seen that the C 1 , C 3 , C 4 , C 5 , C 6 and C 7 curves monotonically increase while C 2 monotonically decreases. As these scores monotonically vary, they do not allow for the identification of the optimal number of features.

Figure 2 : 7 mAlgorithm 1 2 . 1 ε

 27121 Figure 2: State-of-the-art constraint scores with respect to the number of selected features m for the Dermatology database: (a) C 1 m , (b) C 2 m , (c) C 3 m , (d) C 4 m , (e) C 5 m , (f) C 6 m , (g) C 7 m

Figure 3 Figure 3 :

 33 Figure3illustrates the variation of the proposed constraint scores ε S and ε SS with respect to the number of features m in the Dermatology database, from which p = 3 prototypes have been selected for each class. It should also be noted that as for all state-of-the-art scores C c (c = 1, . . . , 7), ε S is normalized by the total number (k • p) 2 of constraints, and ε SS is normalized by the total number n 2 of cells in the similarity matrix W . The curves of ε S and ε SS are both pseudo-convex . The proposed supervised and semi-supervised constraint scores and their minima can thus help used to identify the optimal number of features.It should be emphasized that although this feature selection procedure eval-

Figure 4 :

 4 Figure 4: Accuracy versus number of selected features m by the supervised constraint scores on eight six benchmark databases

Figure 5 :

 5 Figure 5: Accuracy versus number of selected features m by supervised feature selection methods on six high dimensionality databases

Figure 6 :

 6 Figure 6: Accuracy versus number of selected m by the semi-supervised constraint scores for eight benchmark databases

  (bagging) and a random selection of features (random subspaces) strategy for generating different data views. The constraint score C 6 is then used to measure feature relevance on each data replicate, and the score average of all features across all ensemble components is used to rank all features. In SCGS the semi-supervised feature selection is formulated as a optimization problem where the objective function contains three terms. The first term is the empirical loss on labeled data sample, the second term is a discriminant regularization term focusing on the local structure of data reflected by pairwise constraints, and the last term is an unsupervised estimation on intrinsic geometry distribution of the whole training data. The wrapper-type supervised feature selection (SFS-1NN) is used as baseline. For the sake of fairness, the selection and evaluation are performed in the same conditions as for the semisupervised constraint scores. For the unsupervised methods, only the unlabeled training data samples are used for the feature selection and are then labeled by the constrained spectral clustering before being used as prototypes by the nearest neighbor classifier.Figure 7 displays the average accuracies over 100 runs achieved by the above feature selection approaches on high dimensionality datasets compared with the r% most relevant features, when the number of prototypes p is set to 3. The mean and standard deviation of accuracies, computed

Figure 7 :

 7 Figure 7: Accuracy versus number of selected features m by semi-supervised and unsupervised feature selection methods on six high dimensionality databases

Figure 8 :

 8 Figure 8: Influence of the scaling parameter σ versus number of selected features m by the proposed constraint scores ε S and ε SS . Top row: results obtained by using ε S . Bottom row: results obtained by using ε SS

5. 4 . 4 .

 44 Comparison between supervised constraint score ε S and semi-supervised constraint score ε SS In this section, we compare the accuracy provided by the supervised constraint score ε S and the semi-supervised constraint score ε SS . For this purpose, the accuracy of the test dataset is obtained using the same nearest neighbor classifier based only on the k prototype subsets X l (|X l | = p). Figure 10 illustrates (a) ε S (b) ε SS

Figure 9 :

 9 Figure 9: Accuracy-scores ε S and ε SS versus number of selected features for the Dermatology database

Figure 10 :

 10 Figure 10: Comparison of accuracy rates obtained by ε SS and ε S scores for eight benchmark databases using nearest neighbor classifier

Table 1 :

 1 Description of benchmark databases used in experiments

	Database	#Features #Samples #Training/#Test #Classes
	WBCD	09	683	410/273	2
	Image Segmentation	19	2310	1155/1155	7
	WDBC	30	569	376/193	2
	Ionosphere	34	351	176/175	2
	Dermatology	34	366	183/183	6
	Libras Movement	90	360	180/180	15
	Multiple Feature	649	2000	1000/1000	10
	CNAE-9	856	1080	540/540	9
	Yale	1024	165	90/75	15
	ORL	1024	400	200/200	40
	Pie10P	2420	210	110/100	10
	ALLAML	7129	72	37/35	2
	C

2

, then with that achieved by supervised feature selection methods including

Table 2 :

 2 Rank sum of different supervised constraint scores for different numbers of prototypes

	Database	m	p	ε S	C 1	C 2
			2	148	171	161
	WBCD	5	3 * 160 * 169 * 160 *
			4	172	167	158
			2	196	183 192
	Image Segmentation	11	3 * 170 * 209 * 184 *
			4	184	202	182
			2	171	188	171
	WDBC	12	3 4	183	176 194
			2	185	192	185
	Ionosphere	13	3 4	179	192	203
			2	170	194	223
	Dermatology	15	3 4	160	202	223
			2	127	233	232
	Libras Movement	45	3 4	125	223	234
			2	156	216	226
	Multiple Features	100	3 4	116	238	245
			2	141	184	248
	ORL	300	3

* 186 * 194 * 176 * * 184 * 194 * 197 * * 164 * 199 * 227 * * 118 * 225 * 251 * * 133 * 218 * 246 * * 151 * 160 * 278 * 4 100 100 100

Table 3 :

 3 Mean and standard deviations of accuracy over the r% selected features

	Database SFS-1NN ReliefF mRMR C 1 C 2 ε S	WBCD 96.75 ± 2.84 ± 3.28 94.99 ± 3.94 95.63 ± 3.19 95.51 ± 3.49 95.53 ± 3.38	Image Segmentation 60.77 ± 6.82 53.91 ± 7.68 55.55 ± 6.97 52.94 ± 7.03 53.80 ± 6.76 55.43 ± 6.58	WDBC 90.88 ± 5.37 87.29 ± 7.49 86.64 ± 7.99 87.41 ± 7.41 87.88 ± 7.35 87.60 ± 7.54	Ionosphere 77.66 ± 6.82 66.95 ± 8.80 66.56 ± 9.59 66.96 ± 8.72 67.02 ± 8.73 68.02 ± 8.44	Dermatology 90.90 ± 2.80 75.10 ± 5.60 83.96 ± 5.50 83.22 ± 5.10 82.47 ± 5.76 84.00 ± 5.29	Libras Movement 38.06 ± 5.01 32.25 ± 4.32 36.17 ± 4.27 33.41 ± 4.43 32.09 ± 4.53 37.32 ± 4.61	Multiple Feature 94.24 ± 1.01 65.99 ± 5.39 74.85 ± 4.55 75.91 ± 5.32 69.89 ± 4.37 76.97 ± 4.31	CNAE-9 67.96 ± 3.73 34.51 ± 6.06 40.38 ± 5.59 42.38 ± 5.72 36.04 ± 5.65 42.61 ± 5.27	Yale 49.57 ± 4.94 41.17 ± 5.73 48.49 ± 5.49 49.05 ± 5.56 31.75 ± 3.79 50.60 ± 5.49	ORL 70.40 ± 3.07 66.77 ± 3.17 72.75 ± 2.92 64.67 ± 3.50 65.19 ± 2.64 73.84 ± 2.44	Pie10P 59.60 ± 11.58 49.63 ± 7.14 44.90 ± 4.55 46.54 ± 5.30 46.45 ± 8.14 51.39 ± 6.50	ALLAML 71,95 ± 7,65 68.51 ± 13.27 64.50 ± 13.10 68.03 ± 13.20 66.51 ± 14.48 67.82 ± 14.48	Average rank 1.25 4.5 4.42 3.84 4.75 2.33

Table 4

 4 lists the means of N M L and CN M L for eight databases over 100 different sets of prototypes when the number of prototypes p ranges from 2 to 29 4. To interpret N M L with set M , Table 4 also displays its size |M | = k•p 2 -k•p.

Table 4 :

 4 N M and CN M L achieved for eight databasesThis table demonstrates that for each database, N M L decreases when the number of prototypes p increases. N M L also varies strongly from one database to another and tends to be high when the number |M | of available pairs is low. Furthermore, CN M L is close to or exceeds 0.7 for all databases except for the Image Segmentation, Libras Movement and Ionosphere databases.

	Database	p |M | N M L CN M L
		2 04	24037 0.8596
	WBCD	3 12	7956.2 0.8647
		4 24	3915.4 0.8799
		2 14	15234 0.4690
	Image Segmentation	3 42	4922.4 0.4966
		4 84	2408.2 0.5126
		2 04	19638 0.7962
	WDBC	3 12	6392	0.8094
		4 24	3140.9 0.8337
		2 04	5450.1 0.5826
	Ionosphere	3 12	1752.6 0.5962
		4 24	879.9	0.6008
		2 12	534.2	0.8050
	Dermatology	3 36	174.4	0.8504
		4 72	86.5	0.8630
		2 30	81.19	0.4305
	Libras Movement	3 90	25.0	0.5077
		4 180 12.0	0.5684
		2 20	5198.6 0.6269
	Multiple Features	3 60	1698.3 0.6939
		4 120 841.3	0.7401
		2 80	11.1	0.6324
	ORL	3 240 3.4	0.7694
		4 480 1.7	0.8780

Table 5 :

 5 Rank sum of different semi-supervised constraint scores for different number of 192 * 304 * 347 * 250 * 251 * 209 * 155 * 282 * 292 * 160 * 164 * 156 * 299 * 351 * 372 * 303 * 306 * 282 *

	prototypes
	Database	m	p	ε SS	C 3	C 4	C 5	C 6	C 7
			2	254	257	374	246	244 251
	WBCD	5	3 4	206	307	332	253	250	221
			2	179	271	284	211	218	179
	Image Segmentation	12	3 4	176	210	242	192	196	176
			2	282	358	359	305	314	324
	WDBC	16	3 4	308	314	341	316	301 306
			2	262	334	332	321	323	379
	Ionosphere	17	3 4	291	332	333	307	303	434
			2	208	316	367	282	292	482
	Dermatology	15	3 4	177	306	374	309	324	466
			2	132	318	331	409	404	432
	Libras Movement	45	3 4	138	312	323	464	483	313
			2	149	486	369	311	311	364
	Multiple Features	100	3 4	108	375	454	382	379	297
			2	154	246	378	508	518	240
	ORL	300	3

* * * * 283 * 326 * 334 * 314 * 307 * 414 * * 166 * 319 * 343 * 314 * 318 * 498 * * 128 * 306 * 326 * 462 * 458 * 358 * * 119 * 422 * 407 * 378 * 375 * 300 * * 174 * 210 * 416 * 567 * 497 * 189 *

Table 6 :

 6 Mean and standard deviations of accuracy over the r% selected features 14.16 62.74 ± 12.48 78.30 ± 12.35 61.55 ± 11.42 58.50 ± 15.19 68.16 ± 11.36 67.35 ± 13.31

	C 7 ε SS	96.20 ± 4.16 96.34 ± 4.00	50.63 ± 6.86 55.26 ± 5.98	82.41 ± 7.07 83.03 ± 7.29	59.74 ± 8.19 63.15 ± 10.53	75.69 ± 6.49 84.98 ± 5.91	31.84 ± 4.74 38.46 ± 4.13	75.66 ± 6.87 84.22 ± 4.41	19.13 ± 2.28 46.60 ± 6.78	48.31 ± 5.56 51.76 ± 4.76	61.09 ± 3.41 75.16 ± 2.33	49.32 ± 8.49 55.56 ± 7.44	3.25 1.67
	Database LS Spec SCGS EnsCLS C 4	WBCD 95.36 ± 4.74 92.58 ± 7.42 92.98 ± 6.66 92.21 ± 8.06 95.14 ± 5.34	Image Segmentation 50.22 ± 6.51 54.39 ± 6.92 46.35 ± 8.29 38.20 ± 5.56 50.48 ± 7.69	WDBC 83.04 ± 7.74 82.05 ± 7.01 81.68 ± 8.60 73.41 ± 8.80 81.41 ± 7.67	Ionosphere 59.86 ± 8.43 64.06 ± 10.52 57.96 ± 8.95 65.29 ± 11.15 61.81 ± 9.19	Dermatology 70.41 ± 6.33 66.20 ± 6.29 71.07 ± 7.97 64.03 ± 5.60 77.84 ± 7.91	Libras Movement 32.07 ± 3.81 29.42 ± 4.01 32.19 ± 4.56 31.25 ± 3.94 32.59 ± 4.51	Multiple Feature 69.55 ± 4.28 70.60 ± 4.06 49.14 ± 4.82 34.07 ± 3.75 71.76 ± 5.74	CNAE-9 11.50 ± 0.40 14.09 ± 1.67 51.71 ± 6.50 11.39 ± 0.30 43.35 ± 7.28	Yale 48.10 ± 3.56 47.82 ± 4.00 30.84 ± 5.02 29.62 ± 4.27 47.4 ± 4.22	ORL 56.86 ± 2.49 55.38 ± 2.50 50.13 ± 3.02 46.86 ± 2.48 56.35 ± 2.95	Pie10P 54.09 ± 8.80 48.26 ± 7.19 61.36 ± 10.84 64.81 ± 7.67 46.90 ± 9.03	ALLAML 63.25 ± Averaged rank 4.08 4.67 4.33 5.75 4.25

Table 7 :

 7 Run-time (in seconds) of feature selection methodsThe choice of feature number to select is open challenge. It is often specified in advance. The state-of-the-art constraint scores and the feature selection methods used in previous experiments are not able to determine the number of selected features. On the contrary, our scores ε S and ε SS offer a possibility

	Database	m	mRMR EnsCLS SFS-1NN	ε S	ε SS
		20	0.135	3.207	318,66	0.491 202.35
	Multiple Features	50	0.329	5.488	981,91	2.142 521.60
		100	0.651	8.277	2752.06	8.416 1111.0
		20	0.241	1.605	85.478	2.783 6.429
	ORL	50	0.614	2.752	248,34	7.925 18.614
		100	1.190	4.117	602,79	20.78 46.412
		20	1.343	0.067	97.033	0.294 3.603
	ALLAML	50	3.417	0.102	250,72	1.045 9.467
		100	7.021	0.146	535.37	3.560 22.452
	5.4.3. Optimal number of selected features

Table 8 :

 8 Average accuracy and optimal number of selected features

		Constraint score ε S	Constraint score ε SS
	Database	Acc(m * )	Acc( m)	Acc(m * )	Acc( m)
	WBCD	96.48 (9)	96.35 (6)	97.41 (6)	96.90 (7)
	Image Segmentation	58.07 (9)	57.70 (18)	58.47 (7)	58.28 (18)
	WDBC	88.54 (14)	87.97 (25)	86.33 (12)	85.24 (25)
	Ionosphere	68.85 (9)	68.73 (12)	65.57 (18)	62.73 (7)
	Dermatology	91.07 (32)	90.82 (25)	92.14 (22)	91.86 (21)
	Libras Movement	42.03 (89)	41.95 (90)	41.73 (41)	39.73 (85)
	Multiple Feature	82.06 (130) 75.27 (24)	87.85 (63)	84.24 (23)
	CNAE-9	45.68 (49)	44.29 (86)	49.43 (65) 46.44 (172)
	Yale	59.24 (205) 49.05 (61)	61.00 (205) 49.19 (51)
	ORL	80.06 (125) 79.91 (115)	81.79 (136) 81.28 (108)
	Pie10P	53.66 (91)	53.26 (67)	65.58 (34)	59.30 (60)
	ALLAML	70.11 (351) 66.29 (29)	71.20 (132) 68.71 (20)

new score can then evaluate the relevance of each pairwise constraint provided by the user.


Weed detection from multispectral images acquired under uncontrolled illumination conditions

Keywords: Multispectral, Weed detection, Reflectance, Band selection, Supervised classification

In order to optimize the application of herbicides in crop fields, selective spraying is one of the promising solutions being investigated during the last decades. Multispectral cameras provide radiance images with a high spectral resolution allowing to analyze vegetation signatures beyond the visible domain in narrow spectral bands. In this study, we address the problem of weed detection using multispectral imaging. A multishot camera is used to acquire radiance images of vegetation in a greenhouse under skylight. From the acquired radiance images, spectral reflectance is estimated in a row-wise manner as an illumination-invariant signature. Reflectance spectra and texture features are extracted and used to feed a classifier in order to conduct supervised weed detection. Furthermore, the most suited spectral bands for the weed/crop separation are identified and used to specify a multispectral filter array (MSFA) for a single-shot (snapshot) camera that can be embed on an unmanned aerial vehicle (UAV) suitable to make instantaneous acquisitions in the field.

MATERIALS AND METHODS

Snapscan camera

The Snapscan 1 is a multispectral camera manufactured by IMEC that embeds a single matrix sensor, covered by a series of narrow stripes of Fabry-Perot integrated filters. It provides K = 141 spectral channels whose central wavelengths range from λ 1 = 475.1 nm to λ 141 = 901.7 nm with a variable center step (from 0.5 nm to 5 nm). Each filter is associated to 5 adjacent rows of 2048 pixels, and samples a band from the visible (VIS) or near infrared (NIR) spectral domain according to its spectral sensitivity function (SSF) T k (λ) with a full width at half maximum between 2 nm and 10 nm.

Multispectral dataset acquisition

An image acquisition campaign was conducted in a greenhouse under skylight by the Chambre d'Agriculture de la Somme in early April 2019. The targeted plants are beet, wheat and bean for crops, and datura, thistle and goose-foot for weeds. The Snapscan camera observes the vegetation from top at a distance of about 1.50 m from the ground (see Fig. 1(a)). Each image contains a portion of a white diffuser (lambertian surface) at the right border and a Color-Checker (see Fig. 1(b)). This campaign provided more than 300 radiance images of 2048 × 2048 pixels ×141 channels of 10-bit depth. The memory size of each image stored in ENVI format is 3.2 GB. For our tests, we use 100 images that characterize vegetation at different growth stages, and containing either mixed vegetation types (crop and weeds) or a single vegetation type (crop or weeds). The spectral reflectance is estimated from the acquired radiance images following a specific method. The reflectance is used as an illumination-invariant feature. Fig. 1(c) shows a classification result obtained using a gradient boosting (LightGBM 2 ) classifier trained with reflectance spectra extracted from training images acquired at different days.

Further author information: (Send correspondence to ) Amziane Anis, E-mail: anis.amziane@univ-lille.fr Macaire Ludovic, E-mail: ludovic.macaire@univ-lille.fr

  Acquisition setup (b) RGB rendering of a multispectral image Beet (Acc = 91%) Thistle (Acc = 98%) (c) Classification result. Beet is shown in green and Thistle (weed) is shown as red.

Figure 1 .

 1 Figure 1. Outdoor multispectral image acquisition and classification.