Cheng-Lin Liu
email: liucl@nlpr.ia.ac.cn

Handwritten Chinese Character Recognition: Effects of Shape Normalization and Feature Extraction

The field of handwritten Chinese character recognition (HCCR) has seen significant advances in the last two decades, owing to the effectiveness of many techniques, especially those for character shape normalization and feature extraction. This paper reviews the major methods of normalization and feature extraction, and evaluates their performance experimentally. The normalization methods include linear normalization, nonlinear normalization (NLN) based on line density equalization, moment normalization (MN), bi-moment normalization (BMN), modified centroid-boundary alignment (MCBA), and their pseudo-two-dimensional (pseudo 2D) extensions. As to feature extraction, we focus on some effective variations of direction features: chaincode feature, normalization-cooperated chaincode feature (NCCF), and gradient feature. We have compared the normalization methods previously, but in this study, will compare them with better implementation of features. As results, the current methods perform superiorly on handprinted characters, but are insufficient for unconstrained handwriting.

Introduction

Since the first work of printed Chinese character recognition (PCCR) was published in 1966 [START_REF] Casey | Recognition of printed Chinese Characters[END_REF], many research efforts have been contributed, to both printed and handwritten Chinese character recognition (HCCR). Research on online HCCR was started as early as PCCR [START_REF] Stalling | Approaches to Chinese character recognition[END_REF], whereas offline HCCR was started in late 1970s, and has attracted high attention from the 1980s [START_REF] Mori | Research on machine recognition of handprinted characters[END_REF]. Since then, many effective methods have been proposed to solve this problem, and the recognition performance has advanced significantly [START_REF] Hildebrandt | Optical recognition of Chinese characters: advances since[END_REF][START_REF] Umeda | Advances in recognition methods for handwritten Kanji characters[END_REF]. This paper is mainly concerned with offline HCCR, but most methods of offline recognition are applicable to online recognition as well [START_REF] Liu | Online recognition of Chinese characters: the state-ofthe-art[END_REF].

The approaches of HCCR can be roughly grouped into two categories: feature matching (statistical classification) and structure analysis. Based on feature vector representation of character patterns, feature matching approaches usually computed a simple distance measure (correlation matching), say, Euclidean or city block distance, between the test pattern and class prototypes. Currently, sophisticated classification techniques [START_REF] Fukunaga | Introduction to Statistical Pattern Recognition[END_REF][START_REF] Duda | Pattern Classification[END_REF][START_REF] Jain | Statistical pattern recognition: a review[END_REF], including parametric and non-parametric statistical classifiers, neural networks, support vector machines (SVMs), etc., can yield higher recognition accuracies. Nevertheless, the selection and extraction of features remains an important issue. Structure analysis is an inverse process of character generation: to extract the constituent strokes and compute a structural distance measure between the test pattern and class models. Due to its resembling of human cognition and the potential of absorbing large deformation, this approach was pursued intensively in the 1980s and is still advancing [START_REF] Kim | Statistical character structure modeling and its application to handwritten Chinese character recognition[END_REF]. However, due to the difficulty of stroke extraction and structural model building, it is not widely followed.

Statistical approaches have achieved great success in handprinted character recognition and are well commercialized due to some factors. First, feature extraction based on template matching and classification based on vector computation are easy to implement and computationally efficient. Second, effective shape normalization and feature extraction techniques, which improve the separability of patterns of different classes in feature space, have been proposed. Third, current machine learning methods enable classifier training with large set of samples for better discriminating shapes of different classes.

The methodology of Chinese character recognition has been largely affected by some important techniques: blurring [START_REF] Iijima | A theoretical study of the pattern identification by matching method[END_REF], directional pattern match-ing [START_REF] Yasuda | An improvement of correlation method for character recognition[END_REF][START_REF] Yamashita | Classification of handprinted Kanji characters by the structured segment matching method[END_REF][START_REF] Fujisawa | Directional pattern matching for character recognition revisited[END_REF], nonlinear normalization [START_REF] Tsukumo | Classification of handprinted Chinese characters using nonlinear normalization and correlation methods[END_REF][START_REF] Yamada | A nonlinear normalization method for hanprinted Kanji character recognition-line density equalization[END_REF], modified quadratic discriminant function (MQDF) [START_REF] Kimura | Modified quadratic discriminant functions and the application to Chinese character recognition[END_REF], etc. These techniques and their variations or improved versions, are still widely followed and adopted in most recognition systems. Blurring is actually a low-pass spatial filtering operation. It was proposed in the 1960s from the viewpoint of human vision, and is effective to blur the stroke displacement of characters of same class. Directional pattern matching, motivated from local receptive fields in vision, is the predecessor of current direction histogram features. Nonlinear normalization, which regulates stroke positions as well as image size, significantly outperforms the conventional linear normalization (resizing only). The MQDF is a nonlinear classifier suitable for high-dimensional features and large number of classes. Its variations include the pseudo Bayes classifier [START_REF] Kimura | Improvement of handwritten Japanese character recognition using weighted direction code histogram[END_REF], the modified Mahalanobis distance [START_REF] Kato | A handwritten character recognition system using modified Mahalanobis distance[END_REF], etc. This paper reviews the major normalization and feature extraction methods and evaluate their performance in offline HCCR on large databases. The normalization methods include linear normalization (LN), nonlinear normalization (NLN) based on line density equalization [START_REF] Tsukumo | Classification of handprinted Chinese characters using nonlinear normalization and correlation methods[END_REF][START_REF] Yamada | A nonlinear normalization method for hanprinted Kanji character recognition-line density equalization[END_REF], moment normalization (MN) [START_REF] Casey | Moment normalization of handprinted character[END_REF], bi-moment normalization (BMN) [START_REF] Liu | Handwritten Chinese character recognition: alternatives to nonlinear normalization[END_REF], modified centroid-boundary alignment (MCBA) [START_REF] Liu | Global shape normalization for handwritten Chinese character recognition: a new method[END_REF], as well as the pseudo-two-dimensional (pseudo 2D) extensions of them [START_REF] Horiuchi | Two-dimensional extension of nonlinear normalization method using line density for character recognition[END_REF][START_REF] Liu | Pseudo twodimensional shape normalization methods for handwritten Chinese character recognition[END_REF]. These methods have been evaluated previously [START_REF] Liu | Pseudo twodimensional shape normalization methods for handwritten Chinese character recognition[END_REF], but in this study, they will be evaluated with better implementation of features.

Though many features have been proposed for character recognition, we focus on the class of direction histogram features, including chaincode direction feature, normalization-cooperated chaincode feature (NCCF) [START_REF] Hamanaka | Normalization-cooperated feature extraction method for handprinted Kanji character recognition[END_REF], and gradient direction feature. These features have yielded superior performance due to the sensitivity to stroke-direction variance and the insensitivity to stroke-width variance. The gradient direction feature was not paid enough attention until the success of gradient vector decomposition [START_REF] Liu | Handwritten digit recognition: benchmarking of state-of-the-art techniques[END_REF], following a decomposition scheme previously proposed in online character recognition [START_REF] Kawamura | On-line recognition of freely handwritten Japanese characters using directional feature densities[END_REF]. Alternatively, the direction of gradient was quantized into a number of angular regions [START_REF] Srikantan | Gradient-based contour encoder for character recognition[END_REF]. By NCCF, the chaincode direction is taken from the original image instead of the normalized image, but the directional elements are displaced in normalized planes according to normalized coordinates. An improved version of NCCF maps chaincodes into continuous line segments in normalized planes [START_REF] Liu | Handwritten digit recognition: investigation of normalization and feature extraction techniques[END_REF].

In the history, some extensions of direction feature, like the peripheral direction contributivity (PDC) [START_REF] Hagita | Handprinted Chinese characters recognition by peripheral direction contributivity feature[END_REF] and the reciprocal feature field [START_REF] Yasuda | An improved correlation method for handprinted Chinese character recognition in a reciprocal feature field[END_REF], have reported higher accuracy in HCCR when simple distance metric was used. These features, with very high dimensionality (over 1,000), are actually highly redundant. As background features, they are sensitive to noise and connecting strokes. Extending the line element of direction feature to higher-order feature detectors (e.g., [START_REF] Teow | Robust vision-based features and classification schemes for off-line handwritten digit recognition[END_REF][START_REF] Shi | Handwritten numeral recognition using gradient and curvature of gray scale image[END_REF]) helps discriminate similar characters, but the dimensionality also increases rapidly. The Gabor filter, also motivated from vision research, promises feature extraction in character recognition [START_REF] Wang | Gabor filterbase feature extraction for character recognition[END_REF], but is computationally expensive compared to chaincode and gradient features, and at best, perform comparably with the gradient feature [START_REF] Liu | Gabor feature extraction for character recognition: comparison with gradient feature[END_REF].

We evaluate the character shape normalization and direction feature extraction methods on two databases of handwritten characters, ETL9B (Electrotechnical Laboratory, Japan) and CASIA (Institute of Automation, Chinese Academy of Sciences), with 3,036 classes and 3,755 classes, respectively. Two common classifiers, minimum distance classifier and modified quadratic discriminant function (MQDF), are used to evaluate the recognition accuracies.

The purpose of this study is twofold. First, the comparison of major normalization and feature extraction methods can provide guidelines for selecting methods in system development. Second, the results show to what degree of performance the state-ofthe-art methods can achieve. We will show in experiments that the current methods can recognize handprinted characters accurately but perform inferiorly on unconstrained handwriting.

In the rest of this paper, we review major normalization methods in Section 2 and direction feature extraction methods in Section 3. Experimental results are presented in Section 4, and finally, concluding remarks are offered in Section 5.

Shape Normalization

Normalization is to regulate the size, position, and shape of character images, so as to reduce the shape variation between the images of same class. Denote the input image and the normalized image by f (x, y) and g(x ′ , y ′ ), respectively, normalization is implemented by coordinate mapping

x ′ = x ′ (x, y), y ′ = y ′ (x, y). (1)
Most normalization methods use 1D coordinate mapping:

x ′ = x ′ (x), y ′ = y ′ (x). (2)
Under 1D normalization, the pixels at the same row/column in the input image are mapped to the same row/column in the normalized image, and hence, the shape restoration capability is limited. Given coordinate mapping functions (1) or ( 2), the normalized image g(x ′ , y ′ ) is generated by pixel value and coordinate interpolation. In our implementation of 1D normalization, we map the coordinates forwardly from (binary) input image to normalized image, and use coordinate interpolation to generate the binary normalized image. For generating gray-scale normalized image, each pixel is viewed as a square of unit area. By coordinate mapping, the unit square of input image is mapped to a rectangle in the normalized plane, and each pixel (unit square) overlapping with the mapped rectangle is assigned a gray level proportional to the overlapping area [START_REF] Liu | Handwritten digit recognition: investigation of normalization and feature extraction techniques[END_REF].

In the case of 2D normalization, the mapped shape of a unit square onto the normalized plane is a quadrilateral [START_REF] Liu | Pseudo twodimensional shape normalization methods for handwritten Chinese character recognition[END_REF]. To compute the overlapping areas of this quadrilateral with the pixels (unit squares) in the normalized plane, the quadrilateral is decomposed into trapezoids such that each trapezoid is within a row of unit squares. Each within-row trapezoid is further decomposed into trapezoids within a unit square. After generating the normalized grayscale image, the binary normalized image is obtained by thresholding the gray-scale image (fixed threshold 0.5).

In our experiments, the normalized image plane is set to a square of edge length L, which is not necessarily fully occupied. To alleviate the distortion of elongated characters, we partially preserve the aspect ratio of the input image. By aspect ratio adaptive normalization (ARAN) [START_REF] Liu | Handwritten digit recognition: investigation of normalization and feature extraction techniques[END_REF][START_REF] Liu | Aspect ratio adaptive normalization for handwritten character recognition[END_REF], the aspect ratio R 2 of normalized image is a continuous function of the aspect ratio R 1 of input image:

R 2 = sin( π 2 R 1 ). ( 3
)
R 1 is calculated by

R 1 = W 1 /H 1 , if W 1 < H 1 H 1 /W 1 , otherwise (4)
where W 1 and H 1 are the width and height of the input image. The width W 2 and height H 2 of the normalized image are similarly related by the aspect ratio R 2 . If the input image is vertically elongated, then in the normalized plane, the vertical dimension is filled (height L) and the horizontal dimension is centered and scaled according to the aspect ratio; otherwise, the horizontal dimension is filled (width L) and the vertical dimension is centered and scaled. ARAN is depicted in Fig. 1.

Figure 1: Aspect ratio adaptive normalization (ARAN). Rectangle with thick line: occupied area of normalized image.

The normalization methods depend on the coordinate mapping functions, defined by the 1D and pseudo 2D normalization methods as follows.

1D normalization methods

Given the sizes of input and normalized images, the coordinate mapping functions of linear normalization (LN) are simply

x ′ = W2 W1 x, y ′ = H2 H1 y. (5)
Both the linear normalization and line-density-based nonlinear normalization (NLN) methods align the physical boundaries (ends of stroke projections) of input image to the boundaries of normalized image. The coordinate mapping of NLN is obtained by accumulating the normalized line density projections (line density equalization):

x ′ = W 2 x u=0 h x (u), y ′ = H 2 y v=0 h y (v), (6)
where h x (x) and h y (y) are the normalized line density histograms of x axis and y axis, respectively, which are obtained by normalizing the projections of local line densities into unity sum:

       h x (x) = px(x) x px(x) = y dx(x,y)
where p x (x) and p y (y) are the line density projections onto x axis and y axis, respectively, and d x (x, y) and d y (x, y) are local line density functions.

By Tsukumo and Tanaka [START_REF] Tsukumo | Classification of handprinted Chinese characters using nonlinear normalization and correlation methods[END_REF], the local line densities d x and d y are taken as the reciprocal of horizontal/vertical run-length in background area, or a small constant in stroke area. While by Yamada et al. [START_REF] Yamada | A nonlinear normalization method for hanprinted Kanji character recognition-line density equalization[END_REF], d x and d y are calculated considering both background run-length and stroke run-length, and are unified to render d x (x, y) = d y (x, y). The two methods perform comparably but the one of Tsukumo and Tanaka is computationally simpler [START_REF] Umeda | Advances in recognition methods for handwritten Kanji characters[END_REF]. By adjusting the density functions of marginal and stroke areas empirically in Tsukumo and Tanaka's method, we have achieved better performance than the method of Yamada et al. This improved version of NLN is taken in our experiments.

The 1D moment normalization (MN) method (a simplified version of Casey's method [START_REF] Casey | Moment normalization of handprinted character[END_REF]) aligns the centroid of input image (x c , y c ) to the geometric cen-

ter of normalized image (x ′ c , y ′ c ) = (W 2 /2, H 2 /2
), and re-bound the input image according to secondorder 1D moments. Let the second-order moments be µ 20 and µ 02 , the width and height of input image are re-set to δ x = 4

√ µ 20 and δ y = 4 √ µ 02 , respectively. The coordinate mapping functions are then given by

x ′ = W2 δx (x -x c ) + x ′ c , y ′ = H2 δy (y -y c ) + y ′ c . (8)
The bi-moment normalization (BMN) method [START_REF] Liu | Handwritten Chinese character recognition: alternatives to nonlinear normalization[END_REF] aligns the centroid of input image as moment normalization does, but the width and height are treated asymmetric with respect to the centroid. To do this, the second-order moments are split into two parts by the centroid: µ - x , µ + x , µ - y , and µ + y . The boundaries of input image are re-set to [x c -

2 µ - x , x c +2 µ + x ] and [y c -2 µ - y , y c +2 µ + y ]. For the x axis, a quadratic function u(x) = ax 2 + bx + c is used to align three points (x c -2 µ - x , x c , x c + 2 µ +
x ) to normalized coordinates (0, 0.5, 1), and similarly, a quadratic function v(y) is used for the y axis. Finally, the coordinate functions are

x ′ = W 2 u(x), y ′ = H 2 v(y). ( 9
)
The quadratic functions can also be used to align the physical boundaries and centroid, i.e., map (0, x c , W 1 ) and (0, y c , H 1 ) to (0, 0.5, 1). We call this method centroid-boundary alignment (CBA). A modified CBA (MCBA) method [START_REF] Liu | Global shape normalization for handwritten Chinese character recognition: a new method[END_REF] also adjusts the stroke density in central area by combining a sine function with the quadratic functions:

x ′ = W 2 [u(x) + η x sin(2πu(x))], y ′ = H 2 [v(y) + η y sin(2πv(y))]. ( 10
)
The amplitudes of sine waves, η x and η y , are estimated from the extent of the central area, which is defined by the centroids of partial images divided by the global centroid.

Pseudo 2D normalization methods

Horiuchi et al. proposed a pseudo 2D nonlinear normalization (P2DNLN) method by equalizing the line density functions of each row/column instead of the line density projections [START_REF] Horiuchi | Two-dimensional extension of nonlinear normalization method using line density for character recognition[END_REF]. To control the degree of shape deformation, they blurred the line density functions such that the equalization of each row/column is dependent on its neighboring rows/columns. Though this method promises recognition, it is computationally expensive due to the row/column-wise line density blurring.

An efficient pseudo 2D normalization approach, called line density projection interpolation (LDPI), was proposed recently [START_REF] Liu | Pseudo twodimensional shape normalization methods for handwritten Chinese character recognition[END_REF]. Instead of line density blurring and row/column-wise equalization, LDPI partitions the 2D line density map into soft strips. 1D coordinate functions are computed from the density projection of each strip and combined to a 2D function. Specifically, let the width and height of the input image be W 1 and H 1 , the centroid be (x c , y c ), we partition the horizontal line density map d x (x, y) into three horizontal strips:

d i x (x, y) = w i (y)d x (x, y), i = 1, 2, 3. (11)
w i (y) (i = 1, 2, 3) are piecewise linear functions:

       w 1 (y) = w 0 yc-y yc , y < y c , w 2 (y) = 1 -w 1 (y), y < y c , w 2 (y) = 1 -w 3 (y), y ≥ y c , w 3 (y) = w 0 y-yc H1-yc , y ≥ y c , (12)
where w 0 controls the weight of the upper/lower part of line density map. A small value of w 0 renders the interpolated 2D coordinate function close to that of 1D normalization, while a large one may yield excessive deformation. The weight functions with w 0 = 1 are depicted in Fig. 2. The horizontal density functions of three strips are projected onto the x axis:

p i x (x) = y d i x (x, y), i = 1, 2, 3. (13)
The projections are then normalized to unity sum and accumulated to give 1D coordinate functions x ′i (x), i = 1, 2, 3, which are combined to 2D coordinate function by interpolation:

x ′ (x, y) = w 1 (y)x ′1 (x) + w 2 (y)x ′2 (x), y < y c , w 3 (y)x ′3 (x) + w 2 (y)x ′2 (x), y ≥ y c . (14)
Similar to the partitioning of horizontal density, the vertical density map d y (x, y) is partitioned into three vertical strips using weight functions in x axis. The partitioned density functions d i y (x, y), i = 1, 2, 3, are similarly equalized and interpolated to generate the 2D coordinate function y ′ (x, y).

The strategy of LDPI is applied to extend other 1D normalization methods: MN, BMN, and MCBA. The extended versions are called pseudo 2D MN (P2DMN), pseudo 2D BMN (P2DBMN), and pseudo 2D CBA (P2DCBA), respectively. These methods do not rely on the computation of local line density map. Instead, they are directly based on the pixel intensity of character image. As the soft partitioning of line density map in LDPI, the input character image f (x, y) is softly partitioned into three horizontal strips f i x (x, y), i = 1, 2, 3, and three vertical strips f i y (x, y), i = 1, 2, 3. The horizontal strips are projected onto the x axis:

p i x (x) = y f i x (x, y), i = 1, 2, 3. ( 15
)
For P2DMN, the second order moment is computed from the projection of a strip:

µ i 20 = x (x -x i c ) 2 p i x (x) x p i x (x) . ( 16
)
The width of this strip is re-set to δ i x = 4 µ i 20 , which is used to determine the scaling factor of 1D coordinate mapping:

x ′i (x) = W 2 δ i x (x -x i c ) + W 2 2 . ( 17
)
And the 1D coordinate functions of vertical strips are computed from strip projections similarly. For P2DBMN, the second order moment of a horizontal strip is split into two parts at the centroid of this strip: µ i- 20 and µ i+ 20 . The bounds of this strip is re-set to [

x i c -2 µ i- 20 , x i c + 2 µ i+ 20 ]
, which, together with the centroid x i c , are used to estimate the quadratic 1D coordinate mapping function x ′i (x).

The three 1D coordinate functions of vertical strips are computed similarly.

By P2DCBA, from the vertical projection of each horizontal strip f i x (x, y) (i = 1, 2, 3), the centroid coordinate x i c and two partial centroids x i c1 and x i c2 are computed to estimate the parameters of 1D coordinate mapping function x ′i (x). Similarly, 1D coordinate functions y ′i (y) (i = 1, 2, 3) are estimated from the horizontal projections of vertical strips. More details of pseudo 2D normalization can be found in [START_REF] Liu | Pseudo twodimensional shape normalization methods for handwritten Chinese character recognition[END_REF]. Some examples of normalization using nine methods (LN, NLN, MN, BMN, MCBA, LDPI, P2DMN, P2DBMN, and P2DCBA) are shown in Fig. 3. We can see that whereas linear normalization (LN) keeps the original shape (only aspect ratio changed), NLN can effectively equalize the line intervals. The centroid-based normalization methods (MN, BMN, and MCBA) effectively regulate the overall shape (skewness of gravity, balance of inner/outer stroke density). The pseudo 2D methods make the stroke positions more uniform, especially, alleviate the imbalance of width/height and position of character parts.

Direction Feature Extraction

The implementation of direction feature is varying depending on the directional element decomposition, the sampling of feature values, the resolution of direction and feature plane, etc. Considering that the stroke segments of Chinese characters can be approximated into four orientations: horizontal, vertical, left-diagonal and right-diagonal, early works usually decomposed the stroke (or contour) segments into four orientations. Feature extraction from stroke contour has been widely adopted because the contour length is nearly independent on stroke-width variation. The local direction of contour, encoded as a chaincode, actually has eight directions (Fig. 4). Decomposing the contour pixels into eight directions instead of four orientations (a pair of opposite directions merged into one orientation) was shown to significantly improve the recognition accuracy [START_REF] Liu | Handwritten digit recognition: benchmarking of state-of-the-art techniques[END_REF]. This is because separating the two sides of stroke edge can better discriminate parallel strokes. The direction of stroke edge can also be measured by the gradient of image intensity, which applies to gray-scale image as well binary image. The gradient feature has been applied to Chinese character recognition in 8-direction [START_REF] Liu | High accuracy handwritten Chinese character recognition using quadratic classifiers with discriminative feature extraction[END_REF] and 12-direction [START_REF] Liu | Handwritten character recognition using gradient feature and quadratic classifier with multiple discrimination schemes[END_REF].

Direction feature extraction is accomplished in In the following, we first describe the directional decomposition procedures for three types of direction features: chaicode direction feature, normalization-cooperated chaincode feature (NCCF), and gradient direction feature, and then address the sampling of direction planes.

Directional decomposition

Directional decomposition results in a number of direction planes (with same size as the normalize image), f i (x, y), i = 1, . . . , N d . We first describe the procedures for decomposing contour/gradient into eight directions, then extend to 12 directions and 16 directions.

In binary image, a contour pixel is a black point with at least one of its 4-connected neighbors being white. The 8-direction chaincodes of contour pixels can be decided by contour tracing, or more simply, by raster scan [START_REF] Liu | Preprocessing and statistical/structural feature extraction for handwritten numeral recognition[END_REF]. At a black pixel (x, y), denoting the values of 8-connected neighbors in counterclockwise as p k , 0, 1, . . . , 7, with the east neighbor being p 0 . For k = 0, 2, 4, 6, if p k = 0, check p k+1 : if

p k+1 = 1 (chaincode k + 1), f k+1 (x, y) increases by 1; otherwise, if p (k+2)%8 = 1 (chaincode (k + 2)%8), f (k+2)%8 (x, y) increases by 1.
For NCCF, each chaincode in the original image is viewed as a line segment connecting two neighboring pixels, which is mapped to another line segment in a standard direction plane by coordinate mapping. In the direction plane, each pixel (unit square) crossed by the line segment in the main (x or y) direction is given a unit of direction contribution. To exploit the continuous nature of line segment, the strength of line direction falling in a pixel is proportional to the length of line segment falling in the unit square (continuous NCCF [START_REF] Liu | Handwritten digit recognition: investigation of normalization and feature extraction techniques[END_REF]). As in Fig. 5, where a line segment mapped from a chaincode covers four unit squares A, B, C and D. By discrete NCCF, the pixels A and C are assigned a direction unit, whereas by continuous NCCF, all the four pixels are assigned direction strengths proportional to the in-square line segment length. In gradient direction feature extraction, the gradi-ent vector, computed on the normalized image using the Sobel operator, is decomposed into components in eight chaincode directions. The Sobel operator has two masks to compute the gradient components in two axes. The masks are shown in Fig. 6, and the gradient g(x, y) = [g x , g y ] T at location (x, y) is computed by

g x (x, y) = f (x + 1, y -1) + 2f (x + 1, y) + f (x + 1, y + 1) -f (x -1, y -1) -2f (x -1, y) + f (x -1, y + 1), g y (x, y) = f (x -1, y + 1) + 2f (x, y + 1) + f (x + 1, y + 1) -f (x -1, y -1) -2f (x, y -1) + f (x + 1, y -1). ( 18
)
Figure 6: Templates of Sobel gradient operator.

The gradient strength and direction can be computed from the vector [g x , g y ] T . The range of gradient direction can be partitioned into a number (say, 8 or 16) of regions and each region corresponds to a direction plane. More effectively, the gradient vector is decomposed into components in standard directions, following a strategy previously proposed in online character recognition [START_REF] Kawamura | On-line recognition of freely handwritten Japanese characters using directional feature densities[END_REF]. In this scheme, if a gradient direction lies between two standard directions, the vector is decomposed into two components as shown in Fig. 7. The length of each component is assigned to the corresponding direction plane at the pixel (x, y). Fig. 8 shows the direction planes of three decomposition schemes: NCCF, chaincodes of normalization image, and gradient of normalized image. The direction planes are arranged in the order of stroke orientation. We can see that the planes of chaincode directions (third row) are very similar to those of gradient directions. The planes of NCCF, describing the local directions of original image, show some difference. Comparing the original image and the normalized image, the orientation of the righthand stroke, near left-diagonal orientation, deforms to near vertical. Consequently, the direction planes of left-diagonal orientation of NCCF are stronger than those of chaincodes and gradient, while the planes of vertical orientation of NCCF are weaker than those of chaincodes and gradient.

Extention to more directions

The extension of gradient decomposition into more than eight directions is straightforward: simply setting N d standard directions with angle interval 360/N d and typically, with one direction pointing to east, then decompose each gradient vector into two components in standard directions and assign the component lengths to corresponding direction planes. We set N d to 12 and 16.

To decompose contour pixels into 16 directions follows the 16-direction extended chaincodes, which is defined by two consecutive chaincodes. In the weighted direction histogram feature of Kimura et al. [START_REF] Kimura | Improvement of handwritten Japanese character recognition using weighted direction code histogram[END_REF], 16-direction chaincodes are down-sampled by weighted average to form 8-direction planes.

Again, we can determine the 16-direction chaincode of contour pixels by raster scan. At a contour pixel (x, y), when its 4-connected neighbor p k = 0 and the counterclockwise successor p k+1 = 1 or p (k+2)%2 = 1, search the neighbors clockwise from p k until a p j = 1 is found. The two contour pixels, p k+1 or p (k+2)%2 and p j , form a 16-direction chaincode. For example, in Fig. 9, the center pixel has the east neighbor being 0, the north neighbor alone defining the 8-direction chaincode, and defining a 16-direction chaincode together with the southeast neighbor. The 16-direction chaincode can be indexed from a table of correspondence between the code and the difference of coordinates of two pixels forming the code, as shown in Fig. 10. Each contour pixel has a unique 16-direction code. For decomposing contour pixels into 12 directions, the difference of coordinates corresponding to a 16direction chaincode is viewed as a vector (the dashed line in Fig. 9), which is decomposed into components in 12 standard directions as a gradient vector is done. In this sense, the 12-direction code of a contour pixel is not unique. For 12-direction chaincode feature extraction, a contour pixel is assigned to two direction planes, with strength proportional to the component length. For NCCF, the two corresponding direction planes are assigned strengths proportional to the overlapping length of the line segment mapped by coordinate functions, as in Fig. 5.

Blurring and sampling

Each direction plane, with same size as the normalized image, need to be reduced to extract feature values of moderate dimensionality. A simple way is to partition the direction plane into a number of block zones and take the total or average value of each zone as a feature value. Partition of variable-size zones was proposed to overcome the non-uniform distribution of stroke density [START_REF] Yamashita | Classification of handprinted Kanji characters by the structured segment matching method[END_REF], but is not necessary when nonlinear or pseudo 2D normalization is done. Overlapping blocks were taken to alleviate the effect of stroke-position variation on the boundary of blocks [START_REF] Guo | Recognition of handwritten characters using pattern transformation method with cosine function[END_REF], yet a more effective way is to partition the plane into soft zones, which follow the principle of low-pass spatial filtering and sampling [START_REF] Liu | Preprocessing and statistical/structural feature extraction for handwritten numeral recognition[END_REF]. The blurring operation of Iijima [START_REF] Iijima | A theoretical study of the pattern identification by matching method[END_REF] implies spatial filtering without down-sampling.

In implementation of blurring, the impulse response function (IRF) of spatial filter is approximated into a weighted window, which is also called a blurring mask. The IRF is often taken as a Gaussian function:

h(x, y) = 1 2πσ 2 x exp - x 2 + y 2 2σ 2 x . ( 19
)
According to the Sampling Theorem, the variance parameter σ x is related to the sampling frequency (the reciprocal of sampling interval). On truncating the band-width of Gaussian filter, an empirical formula was given in [START_REF] Liu | Preprocessing and statistical/structural feature extraction for handwritten numeral recognition[END_REF]:

σ x = √ 2t x π , (20)
where t x is the sampling interval. At a location (x 0 , y 0 ) of image f (x, y), the convolution gives a sampled feature value

F (x 0 , y 0 ) = x y f (x, y) • h(x -x 0 , y -y 0 ). ( 21
)
Fig. 11 show the blurred images (without downsampling) of the direction planes in Fig. 8. By blurring, the sparse pixels in direction planes merge into strokes or blobs.

For ease of implementation, we partition a direction plane into a mesh of equal-size blocks and set the sampling points to the center of each block. Assume to extract K × K values from a plane, the size of plane is set to Kt x × Kt x . From N d direction Figure 11: Blurred images (not down-sampled) of the direction planes in Fig. 8. planes, the total number of extracted feature values is

N d × K 2 .
The extracted feature values are causal variables. Power transformation can make the density function of causal variables closer to Gaussian [START_REF] Fukunaga | Introduction to Statistical Pattern Recognition[END_REF]. This is helpful to improve the classification performance of statistical classifiers. Power transformation is also called variable transformation [START_REF] Wakabayashi | On the size and variable transformation of feature vector for handwritten character recognition[END_REF] or Box-Cox transformation [START_REF] Heiden | The Box-Cox metric for nearest neighbor classification improvement[END_REF]. We transform each feature value with power 0.5 without attempt to optimize the transformation functions.

Performance Evaluation

We first compare the performance of the normalization and feature extraction methods on two large databases of handprinted characters (constrained writing), and then test the performance on a small set of unconstrained handwritten characters.

We use two classifiers for classification: the Euclidean distance to class mean (minimum distance classifier) and the MQDF [START_REF] Kimura | Modified quadratic discriminant functions and the application to Chinese character recognition[END_REF]. For reducing the classifier complexity and improving classification accuracy, the feature vector is transformed to a lower dimensionality by Fisher linear discriminant analysis (FLDA) [START_REF] Fukunaga | Introduction to Statistical Pattern Recognition[END_REF]. We set the reduced dimensionality to 160 for all feature types.

Denote the d-dimensional feature vector (after dimensionality reduction) by x, the MQDF of class ω i is computed by

g 2 (x, ω i ) = k j=1 1 λij [(x -µ i ) T φ ij ] 2 + 1 δi x -µ i 2 - k j=1 [(x -µ i ) T φ ij ] 2 + k j=1 log λ ij + (d -k) log δ i , (22)
where µ i is the mean vector of class ω i , λ ij and φ ij , j = 1, . . . , d, are the eigenvalues and eigenvectors of the covariance matrix of class ω i . The eigenvalues are sorted in non-ascending order and the eigenvectors are sorted accordingly. k denotes the number of principal axes, and the minor eigenvalues are replaced with a constant δ i . We set a class-independent constant δ i which is proportional to the average feature variance, with the multiplier selected by 5-fold holdout validation on the training data set.

In our experiments, k was set to 40. The classification of MQDF is speeded up by selecting 100 candidate classes using Euclidean distance. The MQDF is then computed on the candidate classes only. Candidate selection is further accelerated by clustering the class means into groups. The input feature vector is first compared to cluster centers, and then compared to the class means contained in a number of nearest clusters. We set the total number of clusters to 220 for the ETL9B database and 250 for the CASIA database.

MQDF is a promising classification method for HCCR. Even higher performance can be achieved by, e.g., discriminative learning of feature transformation and classifier parameters [START_REF] Liu | High accuracy handwritten Chinese character recognition using quadratic classifiers with discriminative feature extraction[END_REF][START_REF] Liu | Handwritten character recognition using gradient feature and quadratic classifier with multiple discrimination schemes[END_REF]. The optimization of classifier, however, is not the concern of this paper.

Performance on handprinted characters

The normalization and feature extraction methods are evaluated on two databases of handprinted characters. The ETL9B database contains the character images of 3,036 classes (71 hiragana, and 2,965 Kanji characters in the JIS level-1 set), 200 samples per class. This database has been widely evaluated by the community [START_REF] Kimura | Improvement of handwritten Japanese character recognition using weighted direction code histogram[END_REF][START_REF] Guo | Recognition of handwritten characters using pattern transformation method with cosine function[END_REF][START_REF] Kato | A handwritten character recognition system using directional element feature and asymmetric Mahalanobis distance[END_REF]. The CASIA database, collected by the Institute of Automation, Chinese Academy of Sciences, in early 1990s, contains the handwritten images of 3,755 Chinese characters (the level-1 set in GB2312-80 standard), 300 samples per class. Some sample images of the CASIA database are shown in Fig. 12.

In the ETL9B database, we use the first 20 and last 20 samples of each class for testing, and the re- We first compare the performance of three direction features with varying direction resolutions with a common normalization method (Tsukumo and Tanaka's nonlinear normalization (NLN) method with our modifications). The direction resolution of features is set to 8, 12, and 16. For each direction resolution, three schemes of sampling mesh are tested. For 8-direction features, the mesh of sampling is set to 7×7 (M1), 8×8 (M2), and 9×9 (M3); for 12-direction, 6 × 6 (M1), 7 × 7 (M2), and 8 × 8 (M3); and for 16-direction, 5 × 5 (M1), 6 × 6 (M2), and 7 × 7 (M3). We control the size of normalized image (direction planes) to be around 64 × 64, and the dimensionality (before reduction) to be less than 800. The settings of sampling mesh are summarized in Table 1. On classifier training and testing using different direction resolutions and sampling schemes, the error rates on the test set of ETL9B database are listed in Table 2, and the error rates on the test set of CA-SIA database are listed in Table 3. In the tables, the chaincode direction feature is denoted by chn, NCCF is denoted by nccf, and the gradient direc- tion feature by grd-g. For chaincode feature extraction, the normalized binary image is smoothed using a connectivity-preserving smoothing algorithm [START_REF] Liu | Preprocessing and statistical/structural feature extraction for handwritten numeral recognition[END_REF]. Gradient feature is extracted from gray-scale normalized image.

We can see that on either database, using either classifier (Euclidean or MQDF), the error rates of 12direction and 16-direction features are mostly lower than those of 8-direction features. This indicates that increasing the resolution of direction decomposition is beneficial. The 16-direction feature, however, does not outperform the 12-direction feature. To select a sampling mesh, let us focus on the results of 12-direction features. We can see that by Euclidean distance classification, M1 (6 × 6) outperforms M2 (7×7) and M3 (8×8), whereas by MQDF, the error rates of M2 and M3 are lower than those of M1. Considering that M2 and M3 perform comparably while M2 has lower complexity, we take the sampling mesh M2 with 12-direction features for fol-lowing experiments. The original dimensionality of direction features is now 12 × 7 × 7 = 588. On fixing the direction resolution (12-direction) and sampling mesh (7 × 7), we combine the three types of direction features with nine normalization methods. The weight parameter w 0 of pseudo 2D normalization was set to 0.75 for good recognition performance [START_REF] Liu | Pseudo twodimensional shape normalization methods for handwritten Chinese character recognition[END_REF]. The error rates on the test sets of two databases are listed in Table 4 and Table 5, respectively. In the tables, the pseudo 2D normalization methods P2DMN, P2DBMN, and P2DCBA are denoted by 2MN, 2BMN, and 2CBA, respectively, for saving space. Comparing the normalization methods, we can see that pseudo 2D methods are superior to 1D ones, and the linear normalization (LN) is inferior to other 1D normalization methods. To compare the three types of features, let us view the error rates of 1D normalization methods and those of pseudo 2D methods separately.

It can be seen in Table 4 and Table 5 that with 1D normalization, the NCCF and the gradient fea-ture perform comparably, and both outperform the chaincode feature. Four normalization methods, namely, NLN, MN, BMN, and MCBA, perform comparably. With pseudo 2D normalization, the NCCF perform best, and the gradient feature outperforms the chaincode feature. Comparing the pseudo 2D normalization methods, LDPI and P2DBMN outperform P2DMN and P2DCBA (especially on the CASIA database). On both two databases, the best performance is yielded by the NCCF with LDPI normalization, and the NCCF with P2DBMN is competing.

The gradient feature performs comparably with NCCF with 1D normalization but inferiorly with pseudo 2D normalization. This is because pseudo 2D normalization, though equalize the stroke density better than 1D normalization, also deforms the stroke directions remarkably. While the gradient feature describes the deformed stroke directions, NCCF takes the stroke directions of original image.

Computation times

To compare the computational complexity of normalization methods, we profile the processing time in two sub-tasks: coordinate mapping and normalized image generation, and the latter is dichotomized into binary image and gray-scale image. Smoothing is involved in binary normalized image, but not for gray-scale image.

On the test samples of CASIA database, we counted the CPU times on Pentium-4-3GHz processor. The average times per sample are shown in Table 6. We can see that the processing time of coordinate mapping varies with the normalization method. The linear normalization (LN) is very fast. NLN is more expensive than other 1D methods, and LDPI is more expensive than other pseudo 2D methods. This is because NLN and LDPI involve line density computation. Nevertheless, all these normalization methods are not computationally expensive. Normalized image generation for pseudo 2D normalization methods is very time consuming because it involves quadrilateral decomposition. The processing time of binary normalized image includes a smoothing time of about 0.3ms.

The CPU time of feature extraction is almost independent of normalization method. It has two parts: directional decomposition and blurring. The average CPU times of three direction features are shown in Table 7. The processing time of blurring depends on the sparsity of direction planes (zero pixels are not considered). The direction planes of chaincode feature are sparse, and those of gra-

Performance on unconstrained handwriting

The error rates of the best methods on handprinted characters, say, 0.64% on ETL9B database and 1.52% on CASIA database, are fairly high. To test the performance on unconstrained handwritten characters, we use the classifiers trained with the training set of CASIA database to classify the samples in a small data set, which contains 20 samples for each of 3,755 classes, written by 20 writers without constraint. Some samples of the unconstrained set are shown in Fig. 13. On the unconstrained set, we evaluate various nor- malization methods with two good features: NCCF and gradient direction feature. The error rates are listed in Table 8. We can see that as for handprinted characters, the pseudo 2D normalization methods yield lower error rates than their 1D counterparts on unconstrained characters as well. However, line-density-based nonlinear normalization methods, NLN and LDPI, are evidently inferior to the centroid-based methods. The P2DBMN method performs best, and the P2DCBA method is competitive. Comparing the two types of direction features, the gradient features outperforms NCCF with 1D normalization methods, but with pseudo 2D normalization, NCCF is superior. Overall, the error rates on unconstrained characters are still very high. Fig. 14 shows some misclassified samples by MQDF with the best normalization-feature combination (P2DBMN and NCCF). Most of the misclassified characters are similar in shape to the assigned class: some are inherently similar (top three rows in Fig. 14), and some others are similar due to cursive writing (fourth and fifth rows). Yet the characters of the bottom row are not highly similar to their assigned classes.

Concluding Remarks

We compared various shape normalization and selected feature extraction methods in offline handwritten Chinese character recognition. For direction feature extraction, our results show that 12-direction and 16-direction features outperform 8-direction feature, and the 12-direction feature has better tradeoff between accuracy and complexity. The comparison of normalization methods shows that pseudo 2D normalization methods outperform their 1D counterparts. On handprinted characters, a line-densitybased method LDPI and a centroid-based method P2DBMN perform best. Comparing three types of direction features, the NCCF and the gradient feature outperform the chaincode feature, compete when 1D normalization is used, but with pseudo 2D normalization, the NCCF outperforms the gradient feature. Overall, the best normalization-feature combinations are LDPI and P2DBMN with NCCF.

We also tested the performance on a small set of unconstrained handwritten characters, with classifiers trained with handprinted samples. The error rates on unconstrained characters are very high, but the comparison of normalization and feature extraction methods reveals new insights. Though pseudo 2D normalization methods again outperform their 1D counterparts, the line-density-based LDPI method is inferior to the centroid-based methods, and the P2DBMN method performs best. The best normalization-feature combination on unconstrained characters is P2DBMN with NCCF.

Training classifiers with unconstrained handwritten samples will be able to improve the accuracy of unconstrained character recognition. To collect large database of unconstrained characters is an urgent task in the near future. To reduce the error rate to a fairly low level (say, 2% on isolated characters), however, simply using the current normalization and feature extraction methods and training current classifiers with larger sample set will not suffice. The methods of shape normalization, feature extraction, and classifier design need to be reconsidered for better recognizing cursively written characters. Training classifiers discriminatively can improve the accuracy of both handprinted and unconstrained character recognition.

Figure 2 :

 2 Figure 2: Weight functions for partitioning line density map into soft strips.

Figure 3 :

 3 Figure 3: Character image normalization by nine methods. The leftmost image is original and the other eight are normalized ones.

Figure 4 :

 4 Figure 4: Eight directions of chaincodes.

Figure 5 :

 5 Figure 5: NCCF on continuous direction plane.

Figure 7 :

 7 Figure 7: Decomposition of gradient vector.

Figure 9 :

 9 Figure 9: 16-direction chaincode formed from two 8-direction chaincodes.

Figure 8 :

 8 Figure 8: Original image and normalized image (top row), 8-direction planes of NCCF (second row), chaincode planes of normalized image (third row), and gradient direction planes (bottom row).

Figure 10 :

 10 Figure 10: Difference of coordinates and the corresponding 16-direction chaincodes.

Figure 12 :

 12 Figure 12: Test samples in the CASIA database.

Figure 13 :

 13 Figure 13: Samples of unconstrained handwritten characters.

Figure 14 :

 14 Figure 14: Misclassified samples (each shown original image and normalized image) of unconstrained handwritten characters.

Table 1 :

 1 Settings of sampling mesh for 8-direction, 12-direction, and 16-direction features.

	Mesh	M1	M2	M3
		zones dim zones dim zones dim
	8-dir 12-dir 6 × 6 432 7 × 7 588 8 × 8 768 7 × 7 392 8 × 8 512 9 × 9 648 16-dir 5 × 5 400 6 × 6 576 7 × 7 784

Table 2 :

 2 Error rates (%) of 8-direction, 12-direction, and 16-direction features on ETL9B database.

			Euclidean			MQDF
	chn	M1	M2	M3	M1	M2	M3
	8-dir	2.91 2.94 3.02 1.08 1.05	1.09
	12-dir 2.43 2.56 2.54 1.02 1.00 0.97
	16-dir 2.52 2.40 2.52 1.20 1.00	1.00
	nccf	M1	M2	M3	M1	M2	M3
	8-dir	2.61 2.61 2.71 0.93 0.87	0.89
	12-dir 2.05 2.06 2.13 0.82 0.77 0.78
	16-dir 2.05 2.04 2.11 0.98 0.85	0.79
	grd-g	M1	M2	M3	M1	M2	M3
	8-dir	2.59 2.58 2.66 0.93 0.90	0.89
	12-dir 2.27 2.30 2.31 0.94 0.86	0.86
	16-dir 2.29 2.19 2.25 1.08 0.94 0.85

Table 3 :

 3 Error rates (%) of 8-direction, 12-direction, and 16-direction features on CASIA database.

			Euclidean			MQDF
	chn	M1	M2	M3	M1	M2	M3
	8-dir	5.95 6.07	6.22 2.45 2.37 2.37
	12-dir 5.34 5.44	5.53 2.34 2.26 2.21
	16-dir 5.44 5.28 5.37 2.72 2.35 2.29
	nccf	M1	M2	M3	M1	M2	M3
	8-dir	5.31 5.35	5.49 1.94 1.92 2.00
	12-dir 4.44 4.48	4.55 1.86 1.75 1.73
	16-dir 4.53 4.42 4.52 2.17 1.92 1.82
	grd-g	M1	M2	M3	M1	M2	M3
	8-dir	5.31 5.34	5.41 2.05 2.01 1.97
	12-dir 4.94 4.90	4.98 2.09 2.00 1.95
	16-dir 4.98 4.85 4.80 2.42 2.09 1.98

Table 4 :

 4 Error rates (%) of various normalization methods on ETL9B database.

			Euclidean			MQDF
		chn	nccf grd-g	chn	nccf grd-g
	LN	6.36	5.94	5.97	2.38	2.09	2.11
	NLN	2.56	2.06	2.30	1.00	0.77	0.86
	MN	2.35	2.07	2.12	0.95	0.83	0.82
	BMN	2.33	2.04	2.09	0.92	0.81	0.80
	MCBA 2.52	2.19	2.27	1.00	0.84	0.86
	LDPI	2.08	1.65	1.90	0.82 0.64 0.73
	2MN	2.05	1.65	1.84	0.86	0.69	0.74
	2BMN 1.97 1.60 1.78	0.84	0.69	0.73
	2CBA	2.13	1.81	1.93	0.86	0.72	0.77

Table 5 :

 5 Error rates (%) of various normalization methods on CASIA database.

			Euclidean			MQDF
		chn	nccf	grdg	chn	nccf grd-g
	LN	11.38 10.46 10.31 4.11	3.49	3.54
	NLN	5.44	4.48	4.90	2.26	1.75	2.00
	MN	5.61	4.89	4.90	2.50	2.04	2.06
	BMN	5.30	4.54	4.56	2.35	1.93	1.92
	MCBA	5.48	4.73	4.83	2.31	1.87	1.96
	LDPI	4.49	3.69	4.21	1.96 1.52 1.75
	2MN	4.99	3.97	4.33	2.24	1.70	1.91
	2BMN	4.63	3.71	4.07	2.15	1.62	1.76
	2CBA	4.75	3.95	4.25	2.11	1.68	1.78

Table 6 :

 6 Average CPU times (ms) of normalization on CASIA database. Binary normalized image involves smoothing.

		coordinate binary grayscale
	LN	0.002	0.318	0.133
	NLN	0.115	0.331	0.143
	MN	0.017	0.321	0.126
	BMN	0.024	0.332	0.135
	MCBA	0.032	0.336	0.137
	LDPI	0.266	1.512	1.282
	P2DMN	0.143	1.514	1.236
	P2DBMN	0.147	1.536	1.274
	P2DCBA	0.156	1.542	1.261
	dient feature are densest. The directional decom-
	position of NCCF is most time-consuming because
	it involves line segment decomposition, and gradi-

ent direction decomposition is more expensive than chaincode decomposition. Overall, NCCF is still the most computationally efficient because it saves the time of normalized image generation. The average CPU time of NCCF ranges from 1.21ms to 1.47ms, covering coordinate mapping and feature extraction. For reference, the average classification time of MQDF (with cluster-based candidate selection) for 3,755 classes is 3.63ms, which can be largely reduced by fine implementation, however.

Table 7 :

 7 Average CPU times (ms) of feature extraction on CASIA database.

		direction blurring
	chn	0.121	0.439
	nccf	0.458	0.752
	grd-g	0.329	1.276

Table 8 :

 8 Error rates (%) on unconstrained handwritten Chinese characters.

		Euclidean	MQDF
		nccf	grdg	nccf	grd-g
	LN	36.04	35.44	21.80	21.46
	NLN	26.61	27.04	16.23	16.96
	MN	25.63	25.18	16.40	16.09
	BMN	25.15	24.53	15.91	15.64
	MCBA	25.53	25.06	15.76	15.66
	LDPI	24.93	25.60	15.12	15.73
	P2DMN	23.48	24.38	14.84	15.37
	P2DBMN 22.75 23.35 14.27 14.64
	P2DCBA	23.80	23.78	14.73	14.79

Acknowledgements

This work is supported in part by the Hundred Talents Program of Chinese Academy of Sciences and the Natural Science Foundation of China (grant no.60543004).


Cheng-Lin Liu
email: liucl@nlpr.ia.ac.cn

Polynomial Network Classifier with Discriminative Feature Extraction

The polynomial neural network, or called polynomial network classifier (PNC), is a powerful nonlinear classifier that can separate classes of complicated distributions. A method that expands polynomial terms on principal subspace has yielded superior performance. In this paper, we aim to further improve the performance of the subspace-featurebased PNC. In the framework of discriminative feature extraction (DFE), we adjust the subspace parameters together with the network weights in supervised learning. Under the objective of minimum squared error, the parameters can be efficiently updated by stochastic gradient descent. In experiments on 13 datasets from the UCI Machine Learning Repository, we show that DFE can either improve the classification accuracy or reduce the network complexity. On seven datasets, the accuracy of PNC is competitive with support vector classifiers.

Introduction

Artificial neural networks (ANNs) with supervised learning have shown superior classification performance in many experiments [START_REF] Holmström | Neural and statistical classifiers-taxonomy and two case studies[END_REF]. Frequently used neural classifiers include the multi-layer Perceptron (MLP), radial basis function (RBF) network, polynomial network, etc. The polynomial network is also known as higher-order neural network (HONN), functional link network, polynomial regression [START_REF] Shürmann | Pattern Classification: A Unified View of Statistical and Neural Approaches[END_REF], or generalized linear discriminant function [START_REF] Fukunaga | Introduction to Statistical Pattern Recognition[END_REF]. In this paper, we call this classifier structure as polynomial network classifier (PNC). Since the outputs of PNC are the weighted combinations of higher-order nonlinear functions of input features, it is powerful to separate pattern classes of complicated distributions.

The PNC can be viewed as a single-layer neural network with the input features and their polynomial terms as the network inputs. For d features, the total number of polynomial terms up to r-th order is [START_REF] Shin | Ridge polynomial networks[END_REF]

D = r i=0 d + i -1 i = d + r r . (1)
With large d, the polynomial network will suffer from high computation complexity and will give degraded generalization performance. The complexity can be reduced by either reducing the number of input features or selecting expanded polynomial terms [START_REF] Shürmann | Pattern Classification: A Unified View of Statistical and Neural Approaches[END_REF]. The former way is more computationally efficient and performs fairly well in practice. A PNC with dimensionality reduction by principal component analysis (PCA) has shown superior performance to multilayer neural networks in previous experiments [START_REF] Kreßel | Pattern classification techniques based on function approximation, Handbook of Character Recognition and Document Image Analysis[END_REF][START_REF] Liu | Handwritten digit recognition: benchmarking of state-of-the-art techniques[END_REF].

On the other hand, constrained polynomial structures with moderate complexity have been proposed, like the pi-sigma network (PSN) [START_REF] Shin | The Pi-sigma network: an efficient higher-order neural network for pattern classification and function approximation[END_REF], the ridge polynomial network (RPN) [START_REF] Shin | Ridge polynomial networks[END_REF], and the reduced multivariate polynomial model (RMPM) [START_REF] Toh | Benchmarking a reduced multivariate polynomial pattern classifier[END_REF]. The general HONN is a sigma-pi network in that it combines the products of features. Rather, the output of a PSN is the product of weighted combinations of features. Its number of weights is thus linear with the number of summation units (the order of polynomials). The output of RPN is the summation of pi-sigma units of different orders, and the order can be increased incrementally. The RMPM combines the univariate polynomials, the polynomial of sum of features and its product with the weighted sum of features. These networks actually involve all the polynomial terms of input features up to certain order, but the weights of polynomials are highly constrained. They hence need polynomials of fairly high order (say, 5 or 6) to approximate complicated functions, and cannot guarantee the precision of approximation in difficult cases.

The PNC with full polynomials on reduced features still have higher complexity than the above constrained networks, but usually, a low order (say, 2 or 3) can achieve a reasonable precision of function approximation. The behavior of a lower-order network on feature subspace is easy to explain and to control. Nevertheless, its performance largely depends on the technique of feature selection or dimensionality reduction. Supervised subspace learning methods, like the Fisher linear discriminant analysis (LDA) [START_REF] Fukunaga | Introduction to Statistical Pattern Recognition[END_REF] and heteroscedastic discriminant analysis [START_REF] Brunzell | Feature reduction for classification of multidimensional data[END_REF][START_REF] Loog | Linear dimensionality reduction via a heteroscedastic extension of LDA: the Chernoff criterion[END_REF], may lead to better separability than the unsupervised PCA. These methods, nevertheless, are based on parametric density assumptions and the learning criterion is only loosely connected to classification error.

In this paper, we propose a subspace-feature-based PNC with discriminative feature extraction (DFE). With any classifier structure, DFE optimizes the subspace parameters together with the classifier parameters under a classificationrelated objective on a training sample set [START_REF] Biem | Pattern recognition using discriminative feature extraction[END_REF]. The subspace thus learned is totally classification-oriented and the subspace learning and classifier learning are best fitted. Overfitting can be overcome by adjusting the dimensionality of subspace and the order of classifier. DFE is mostly based on the minimum classification error (MCE) criterion of Juang and Katagiri [START_REF] Juang | Discriminative learning for minimum error classification[END_REF], and has been successfully applied to many pattern recognition problems [START_REF] Wang | Feature extraction and dimensionality reduction algorithms and their applications in vowel recognition[END_REF][START_REF] Yang | Discriminative training approaches to fabric defect classification based on wavelet transform[END_REF]. It has not been combined with polynomial networks, however. Despite that the MCE criterion is applicable to any classifier structures, for neural networks with sigmoid outputs, the minimum squared error (MSE) criterion works well and is easy to optimize by stochastic gradient descent [START_REF] Robbins | A stochastic approximation method[END_REF].

We have evaluated the classification performance of PNC on 13 datasets from the UCI Machine Learning Repository [START_REF]UCI Machine Learning Repository[END_REF]. The results show that compared with the PNC with PCA, DFE either improves the classification accuracy or reduces the network complexity. The complexity of PNC is much lower than support vector classifiers (SVCs) [START_REF] Burges | A tutorial on support vector machines for pattern recognition[END_REF], and on seven of the 13 datasets, the PNC with DFE competes with SVCs in accuracy.

2 Subspace-Feature-Based PNC

We consider second-order (binomial) and third-order polynomial networks, and to save space, we only give the details of binomial networks. The structure and the learning algorithm of third-order networks are similar to binomial ones.

For M -class classification, the PNC has M output units. On a d-dimensional feature vector x = [x 1 , . . . , x d ] T , the output of binomial network for class ω k is computed by

y k (x) = σ d i=1 d j=i w (2) kij x i x j + d i=1 w (1) ki x i + w k0 , (2)
where σ(a) is the sigmoid activation function:

σ(a) = 1 1 + e -a .
In classification, the input pattern (feature vector) is classified to the class of maximum output. The sigmoid function is used in training, and is not necessary in classification. Without the sigmoid function, the network weights can also be estimated by (non-iterative) pseudo inverse [START_REF] Shürmann | Pattern Classification: A Unified View of Statistical and Neural Approaches[END_REF]. Since the sigmoid function makes the network outputs approximate posterior class probabilities, the trained weights with it are more suitable for classification than for regression.

By principal component analysis (PCA), the feature vector is projected onto an m-dimensional principal subspace (m < d):

z = Φ T x = [φ T 1 x, . . . , φ T m x] T = [z 1 , . . . , z m ] T , (3)
where Φ = [φ 1 , . . . , φ m ] is the transformation matrix (subspace basis) composed of the eigenvectors of covariance matrix E[xx T ] corresponding to the m largest eigenvalues. We assume that the origin of the feature space has been shifted to the mean of samples. On the subspace features, the network outputs are computed by

y k (x) = σ m i=1 m j=i w (2)
kij z i z j + m i=1 w (1)
ki z i + w k0 . (4)
On a training set of N samples (x n , c n ), n = 1, . . . , N (c n is the class label of x n ), the connecting weights of PNC are adjusted to minimize the regularized squared error:

E = 1 2N N n=1 M k=1 y k (x n ) -t n k 2 + β w∈W w 2 , ( 5
)
where β is a coefficient of weight decay (excluding the biases); t n k is the target value of class k, with value 1 for the genuine class and 0 otherwise. The weights and biases are initialized to small random values, and by stochastic gradient descent, they are iteratively updated on the training samples until the squared error approaches the minimum. In training, the subspace basis Φ remains unchanged, and the polynomials of projected features can be viewed as the inputs of a single-layer network, for which the training process converges fast.

PNC with Discriminative Feature Extraction

A problem with the subspace-feature-based PNC is that the subspace does not necessarily lead to optimal classification because it is learned independently of the network weights. The subspace learned by PCA does not even consider the class information of training samples. Supervised subspace learning techniques, like LDA and heteroscedastic discriminant analysis, are expected to give better separability than PCA, but do not guarantee the optimality. We aim to learn a better subspace for PNC using discriminative feature extraction (DFE) [START_REF] Biem | Pattern recognition using discriminative feature extraction[END_REF].

By DFE, we adjust not only the network weights in supervised learning, but also the subspace basis simultaneously. Consider that z i = φ T i x, i = 1, . . . , m, let us re-write the network outputs of (4) as

y k (x) = σ m i=1 m j=i w (2) kij φ T i xφ T j x + m i=1 w (1)
ki φ T i x + w k0 = σ(s k (x)), (6)
where s k (x) denotes the weighted sum of output unit k.

In the PNC with DFE, since the projected feature

z i = φ T i x = d j=1 φ ij
x j is a weighted combination of original features and the weights (subspace parameters φ ij , j = 1, . . . , d) are now adjustable, an m-th order polynomial as m i=1 (φ T i x) is actually a pi-sigma unit of the ridge polynomial network (RPN). However, our network has more polynomial terms and needs a lower order than the RPN. Interpreting φ i , i = 1, . . . , m, as subspace basis vectors or feature extractors, a lower-order polynomial network on this feature subspace has decision boundaries of moderate complexity.

The network weights and the subspace basis parameters are adjusted to minimize the regularized square error (5) on a training sample set. The subspace parameters can be initialized to small random values as the network weights. Alternatively, the subspace learned by PCA or LDA is a good start of parameter search. The weights and subspace parameters are then adjusted by stochastic gradient descent on training samples. At time t, the parameters are adjusted on a training sample x by

w (2) kij (t + 1) = w (2) kij (t) -ǫ(t) (y k -t k )y k (1 -y k )z i z j + β N w (2) kij (t) , w (1)
ki (t + 1) = w (1) ki (t) -ǫ(t) (y k -t k )y k (1 -y k )z i + β N w (1) ki (t) , w k0 (t + 1) = w k0 (t) -ǫ(t)(y k -t k )y k (1 -y k ), φ i (t + 1) = φ i (t) -ǫ(t) M k=1 (y k -t k )y k (1 -y k ) ∂s k ∂φi , k = 1, . . . , M, i = 1, . . . , m, j = i, . . . , m, (7)
where ǫ(t) is the learning step, which is set to a small value initially and decreases gradually in the training process. The partial derivative of φ i is specified as

∂s k ∂φ i = 2w (2)
kii z i + j<i w (2) kji z j + j>i w (2) kij z j + w (1) ki x. (8)
In discriminative learning, we keep the unit norm of basis vectors but not the orthogonality. On adjusting the basis vectors on a training sample, each vector is normalized to unit norm ( φ i = 1).

By stochastic gradient descent, the training samples are fed to the PNC for a number (40 or more in our experiments) of cycles. The learning step decreases linearly until it vanishes at the end of training. On every input sample, the network weights and subspace parameters are updated according to [START_REF] Shin | The Pi-sigma network: an efficient higher-order neural network for pattern classification and function approximation[END_REF]. The network weights and the subspace vectors have remarkably different magnitudes of derivatives. To accelerate the convergence of training, they are set two different learning steps, ǫ 1 for weights and ǫ 2 for subspace vectors, and ǫ 1 ≫ ǫ 2 holds.

Another factor affecting training convergence and classification performance is the scale of projected features. We normalize the scale with the square root of the largest eigenvalue λ 1 of E[xx T ] (estimated on training samples and fixed):

z i = φ T i x √ λ 1 . (9)
All the feature vectors are subtracted from the mean of the training samples.

For datasets that have significantly different scales among feature dimensions, it is helpful to uniform the standard deviation of all dimensions of training data (and test data accordingly). This is done before subspace projection.

Experiments

We evaluated the classification performance of subspace-feature-based PNC on 13 datasets from the UCI Machine Learning Repository [START_REF]UCI Machine Learning Repository[END_REF], as summarized in Table 1. We selected the multi-class datasets that have at least 10 features. Some data sets have been partitioned into standard training and test subsets.

For the others, we arrange the samples in random order and evaluate in 5-fold cross-validation. Some datasets have appreciable variability of scale among different dimensions. We normalized them by dividing each dimension with (0.9σ 2 i + 0.1σ 2 0 ) 1/2 , where σ 2 i is the dimension-wise variance and σ 2 0 is the average variance, both estimated on training data.

We compare the PNC-DFE (PNC combined with DFE) with PNC-PCA, oneversus-all support vector classifiers with polynomial and RBF kernels (SVC-poly and SVC-rbf), and the k-nearest neighbor (k-NN) classifier. For the SVC-poly, the feature vectors are uniformly scaled such that the average self-inner product of training vectors is one, and so, the kernel k(x 1 , x 2 ) = (1 + κx 1 • x 2 ) r with κ = 2 i performs fairly well. For the SVC-rbf, the average within-class variance is scaled to one, such that in the kernel function

k(x 1 , x 2 ) = exp(-x1-x2 2 2σ 2
), a parameter value of σ 2 = 0.5 × 2 i performs fairly well. For both, i is an integer selected from -4 to 4.

For the k-NN classifier, SVC-poly, and SVC-rbf, we tried several values of k, polynomial order and κ, or σ 2 such that the classification accuracy on each test set is maximized.

For PNC-PCA and PNC-DFE, we set the number of subspace features m = l • m 1 , l = 1, . . . , 5. m 1 is dependent on the dataset. The selected values of polynomial order and m 1 are listed in the right columns of Table 1. As seen, three datasets (Segment, Pendigit and Letter) are used 2nd-order and the others are used 3rd-order.

The test accuracies (%) of PNC (with full polynomials and dimensionality reduction by PCA and DFE) on the 13 datasets are shown in Table 2. For the "Vowel" dataset, there is no dimensionality reduction when m = 10. For each dataset, the accuracy of full PNC is shown below the title of dataset, and the accuracies of PNC-PCA and PNC-DFE with variable subspace dimensionality are listed in two rows. For the "Isolet" dataset, we do not give the accuracy of full PNC because the number of features is too large.

We can see that on four datasets (Vehicle, Segment, Satimage, and Letter), the full PNC gives the highest accuracy. This can be explained that the four datasets have small number of features and are difficult to classify, so dimensionality reduction by either PCA or DFE cannot improve the classification accuracy. For the other datasets, except for "Soybean-small" and "Isolet", subspace-feature-based PNC performs significantly better than the full PNC.

Comparing the accuracies of PNC-PCA and PNC-DFE, it is evident that except for two datasets (Waveform and Satimage), PNC-DFE mostly give higher accuracies than PNC-PCA, especially on subspaces of lower dimensionality. On The highest accuracies of PNC (full and PNC-PCA), PNC-DFE, SVC-poly, SVC-rbf, and k-NN classifier on the 13 datasets are compared in Table 3. On the "Soybean-small" dataset, all these classifiers achieves perfect classification. Among the other datasets, SVC-poly or SVC-rbf gives the highest accuracies on seven datasets, and PNC or PNC-DFE performs best on five datasets. Expect for four datasets (Soybean-small, Segment, Satimage, and Letter), PNC or PNC performs significantly better than the k-NN classifier. The accuracy of PNC or PNC-DFE is comparable or higher than SVC on seven datasets (Waveform, Wine, Soybean-small, Vehicle, Thyroid, Optdigit, Vowel).

We did not implement the reduced multivariate polynomial model (RMPM) [START_REF] Toh | Benchmarking a reduced multivariate polynomial pattern classifier[END_REF], but results on 10 of our 13 datasets were reported in the literature. Though The complexity of PNC mainly depends on the number of features, and is much lower than SVC and k-NN classifier. The k-NN classifier stores all training samples and compares them with each test pattern. The SVC has a large number of support vectors, ranging from 10% to 70% of all training samples. Due to the limited space, we do not discuss the computational complexity in details.

Conclusion

We proposed to improve the performance of subspace-feature-based polynomial network classifier (PNC) using discriminative feature extraction (DFE), which optimizes the subspace parameters together with the network weights on training samples. Under a regularized squared error criterion, the parameters are efficiently adjusted by stochastic gradient descent. In our experiments on 13 datasets of UCI Machine Learning Repository, DFE mostly improves the accuracy of subspace-feature-based PNC. At moderate complexity, the PNC (full or subspace-feature-based) outperforms the k-NN classifier on nine datasets and competes with support vector classifiers on seven datasets.

Table 1 .

 1 Summary of 12 datasets from UCI Repository. The right two columns shows the selected polynomial order and subspace dimensionality (multiple of m 1).

	Name	#class #feature #train #test Normal. Order m1
	Waveform	3	21	50,000 5-fold No	2	1
	Wine	3	13	178 5-fold Yes	2	2
	Soybean-small	4	35	47 5-fold Yes	2	1
	Vehicle	4	18	846 5-fold Yes	2	2
	Dermatology	6	34	358 5-fold Yes	2	2
	Segment	7	19	2,310 5-fold Yes	3	3
	Thyroid	3	21	3,772 3,428 Yes	2	4
	Satimage	6	36	4,435 2,000	No	2	5
	Optdigit	10	64	3,823 1,797	No	2	10
	Pendigit	10	16	7,494 3,498	No	3	3
	Vowel	11	10	528	462	No	2	2
	Isolet	26	617	6,238 1,559	No	2	25
	Letter	26	16	16,000 4,000	No	3	3

Table 2 .

 2 Test accuracies (%) of PNC (full and PCA) and PNC-DFE on 12 datasets.

	Dataset	PCA=m1 PCA=2m1 PCA=3m1 PCA=4m1 PCA=5m1
	Full PNC DFE=m1 DFE=2m1 DFE=3m1 DFE=4m1 DFE=5m1
	Waveform	63.78	87.02	87.22	87.12	86.98
	84.92	60.64	86.90	86.96	86.72	86.64
	Wine	77.53	90.45	92.13	92.70	92.70
	92.13	79.79	92.13	93.82	93.82	93.82
	Soybean-small 74.47	100	100	100	100
	100	91.49	100	100	100	100
	Vehicle	53.19	67.38	71.75	77.30	78.37
	84.16	71.51	77.42	78.72	79.67	80.02
	Dermatology	77.65	89.94	96.37	96.37	96.37
	96.37	93.30	96.93	96.93	96.65	96.37
	Segment	61.08	84.16	92.21	92.38	92.25
	96.41	92.90	94.33	94.89	95.50	95.80
	Thyroid	92.65	93.49	93.17	93.49	95.51
	94.78	96.06	97.32	97.67	97.72	97.87
	Satimage	86.75	87.10	87.75	88.00	87.95
	88.65	86.45	87.45	87.90	88.15	88.05
	Optdigit	95.72	98.05	98.61	98.61	98.55
	98.50	97.16	98.50	98.72	98.50	98.66
	Pendigit	85.11	95.77	97.68	98.03	98.37
	98.23	89.57	97.17	97.91	98.17	98.37
	Vowel	43.72	50.09	58.01	60.39
	59.52	57.14	60.17	64.72	61.47
	Isolet	93.33	95.19	95.51	96.28	96.28
		95.57	95.96	95.96	96.28	96.09
	Letter	32.35	73.17	85.38	91.47	94.03
	94.70	57.00	80.88	89.47	92.70	94.40
	some datasets (Waveform, Soybean-small, Dermatology, Optdigit, Isolet), the
	PNC-DFE achieves the best or nearly best accuracy on a very low-dimensional
	subspace as m = 2m 1 .

Table 3 .

 3 Highest accuracies of PNC (full and PCA), PNC-DFE, SVCs and k-NN classifier. the datasets were partitioned in different ways, nine of the 10 best accuracies of RMPM (Waveform 83.3%, Soybean-small 95.0%, Vehicle 82.3%, Segment 94.1%, Thyroid 94.0%, Satimage 88.2%, Optdigit 95.3%, Pendigit 95.7%, Letter 74.1%) are lower than our best accuracies of PNC or PNC-DFE.

		PNC PNC-DFE SVC-poly SVC-rbf k-NN
	Waveform 87.22 86.96	87.14	87.08 85.24
	Wine	92.70 93.82	92.13	93.26 87.08
	Soybean-small 100	100	100	100	100
	Vehicle	84.16 80.02	81.56	81.21 71.99
	Dermatology 96.37	96.93	97.77	97.21 96.09
	Segment	96.41	95.80	96.62	96.88 96.71
	Thyroid	95.51 97.87	96.70	95.36 94.28
	Satimage	88.65	88.15	90.70	91.40 90.35
	Optdigit	98.61	98.72	98.66	98.89 98.00
	Pendigit	98.37	98.37	98.77	98.74 97.80
	Vowel	60.39 64.72	56.06	64.50 59.52
	Isolet	96.28	96.28	96.92	96.86 92.69
	Letter	94.70	94.40	96.78	97.65 95.83

Acknowledgements

This work is supported by the Hundred Talents Program of Chinese Academy of Sciences. The author thanks the anonymous reviewers for valuable comments.