IMAGE PROCESSING APPARATUS, IMAGE PROCESSING METHOD, AND COMPUTER PROGRAM PRODUCT

Abstract

An image processing apparatus includes an operating unit that calculates a tangent-line direction of an edge and a normal-line direction of the edge by using a vertical-direction derivative value and a horizontal-direction derivative value of a pixel value of a pixel in an input image; a converting unit that converts local coordinates positioned within a predetermined area with respect to a target pixel in the input image into rotated coordinates, by rotating the local coordinates according to an angle formed by the horizontal direction and the tangent-line direction of the edge; and a fitting unit that performs, at the rotated coordinates, a fitting process that employs a least-squares method by using a curved-plane model expressed with the pixel value of the target pixel in the input image.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is based upon and claims the benefit of priority from the prior Japanese Patent Application No. 2008-222696, filed on Aug. 29, 2008; the entire contents of which are incorporated herein by reference.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The present invention relates to an image processing apparatus, an image processing method, and an image processing computer program product.

2. Description of the Related Art

The specification in the publication of patent application of United States No. 2007/0047838 (reference 1), discloses a technique for obtaining an ideal noise-free image surface by performing a least-squares fitting process while using a parametric curved plane. According to this technique, a method called “kernel regression” is used. With this method, it is possible to extract shapes and edges in an input image with a higher level of precision than when an ordinary linear filtering method is used. However, problems remain where the image becomes too blurry due to a smoothing process when a low-order polynomial curved plane is used, whereas the image is more likely to have noises due to over-fitting when a high-order polynomial curved plane is used.

SUMMARY OF THE INVENTION

According to one aspect of the present invention, an image processing apparatus includes a calculating unit that calculates at least one of a tangent-line direction and a normal-line direction of an edge around a target pixel in an input image by using a direction and a magnitude of a gradient of pixel values of pixels that neighbor the target pixel in the input image; a converting unit that converts coordinates of pixels around the target pixel into a rotated coordinate system by rotating the axes according to an angle between a horizontal direction of the input image and the tangent-line direction of the edge or the normal-line direction of the edge; and a fitting unit that performs a fitting process to fit a local gradient around the target pixel, on the rotated coordinate system, into a curved-plane model by the method of least-squares, and an outputting process to obtain a reconstructed pixel value of the target pixel by using the curved-plane model.

According to another aspect of the present invention, an image processing method includes calculating at least one of a tangent-line direction and a normal-line direction of an edge around a target pixel in an input image by using a direction and a magnitude of a gradient of pixel values of pixels that neighbor the target pixel in the input image; converting coordinates of pixels around the target pixel into a rotated coordinate system by rotating the axes according to an angle between a horizontal direction of the input image and the tangent-line direction of the edge or the normal-line direction of the edge; and performing a fitting process to fit a local gradient around the target pixel, on the rotated coordinate system, into a curved-plane model by the method of least squares, and an outputting process to obtain a reconstructed pixel value of the target pixel by using the curved-plane model.

A computer program product according to still another aspect of the present invention causes a computer to perform the method according to the present invention.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram illustrating an image processing apparatus according to a first embodiment of the present invention;

FIG. 2 is a diagram illustrating an image processing unit according to the first embodiment;

FIG. 3 is a flowchart of operations performed by the image processing unit according to the first embodiment;

FIG. 4 is a drawing for explaining a relationship between a two-dimensional Gaussian function and the tangent-line direction and the normal-line direction of an edge;

FIG. 5 is a drawing for explaining a rotation angle formed by an x-axis of an image and a long-axis direction of a two-dimensional Gaussian function;

FIG. 6 is a flowchart of operations in a curved-plane fitting process performed by the image processing unit according to the first embodiment;

FIG. 7A is a drawing illustrating an original image;

FIG. 7B is a drawing illustrating the shapes of two-dimensional Gaussian functions at points in the original image;

FIGS. 8A and 8B are drawings illustrating results of a fitting process performed with respect to x-y axes according to a conventional kernel regression method;

FIGS. 8C and 8D are drawings illustrating results of a fitting process performed with respect to u-v axes;

FIG. 9 is a drawing for explaining differences between the conventional method and a method according to the first embodiment;

FIG. 10 is a diagram illustrating an image processing unit according to a second embodiment of the present invention;

FIG. 11 is a flowchart of operations performed by the image processing unit according to the second embodiment;

FIG. 12 is a drawing for explaining an image feature classification by Harris;

FIG. 13 is a diagram illustrating an image processing unit according to a third embodiment of the present invention;

FIG. 14 is a flowchart of operations performed by the image processing unit according to the third embodiment;

FIG. 15 is a diagram illustrating an image processing unit according to a fourth embodiment of the present invention;

FIG. 16 is a flowchart of operations performed by the image processing unit according to the fourth embodiment;

FIG. 17 is a diagram illustrating an image processing unit according to a fifth embodiment of the present invention; and

FIG. 18 is a flowchart of operations performed by the image processing unit according to the fifth embodiment.

DETAILED DESCRIPTION OF THE INVENTION

Exemplary embodiments of the present invention will be explained, starting with a first embodiment of the present invention. As shown in FIG. 1, an image processing apparatus 100 includes an input interface unit 101, an image processing unit 102, a storage unit 103, and an output interface unit 104.

The input interface unit 101 is connected to an input device that inputs an image to the image processing apparatus 100. The input interface unit 101 obtains the input image from the input device. The image processing unit 102 performs processes such as a filtering process on the input image. By performing the processes, the image processing unit 102 generates an image from which noises have been eliminated. The storage unit 103 stores therein, for example, the input image that has been obtained by the input interface unit 101 and the image that has been processed by the image processing unit 102. The output interface unit 104 is connected to an output device that outputs images. The output interface unit 104 outputs the images that have been stored in the storage unit 103 to the output device.

One or both of the input device and the output device may be provided on the outside of the image processing apparatus 100 or on the inside of the image processing apparatus 100.

As shown in FIG. 2, an image processing unit 102a that corresponds to the image processing unit 102 shown in FIG. 1 includes an calculating unit 201, a parameter calculator 202, a converting unit 203, and a fitting unit 204.

Next, operations performed by the image processing unit 102a shown in FIG. 2 will be explained with reference to FIG. 3. In the following sections, a position in an input image will be expressed as x (where xεΩ), whereas a set made up of the entirety of positions in the input image will be expressed as Ω (where Ω⊂R²), while a pixel value in the position x in the input image will be expressed as I(x). The pixel value may be, for example, a scalar value (e.g., a luminance value) or a vector value (e.g., a Red/Green/Blue [RGB] color signal). A transposition of a matrix or a vector will be expressed by adding a superscript “T” to the matrix or the vector. Num(N) denotes the number of elements in a set N.

At step S301, for each of the pixels in the input image, the calculating unit 201 calculates a direction and a magnitude of the gradient of a pixel value in a surrounding area of the pixel. According to the first embodiment, an example will be explained in which the calculating unit 201 performs a discretized first-order derivative operation to calculate the direction and the magnitude of the gradient of the pixel value. In the discretized first-order derivative operation, it is possible to use, for example, a Sobel operator. In the case where a Sobel operator is used, the calculating unit 201 performs calculations as shown in Expressions (1) below.

$\begin{array}{cc}{d}_x\left(x\right)=\left[\begin{array}{c}-101\\ -202\\ -101\end{array}\right]*I\left(x\right)d_y\left(x\right)=\left[\begin{array}{ccc}-1& -2& -1\\ 0& 0& 0\\ 1& 2& 1\end{array}\right]*I\left(x\right)& \left(1\right)\end{array}$

In Expressions (1), d_x(x) denotes an x-direction derivative value in the position x, whereas d_y(x) denotes a y-direction derivative value in the position x. The symbol “*” denotes a convolution operation. It is possible to express the value d_x (x) by using Expression (2) shown below.

d_x(x)=−I(x+(−1,−1)^T)+I(x+(1,−1)^T)−2I(x+(−1,0)^T)+2I(x+(1,0)^T)−I(x+(−1,1)^T)+I(x+(1,1)^T)  (2)

It is acceptable to perform a discretized derivative operation that uses any other elements besides the Sobel operator. For example, it is acceptable to use any of the followings: a forward difference, a backward difference, and a central difference.

At step S302, the parameter calculator 202 calculates image feature parameters indicating the direction and the magnitude of an edge in the position x, by using d_x(x) and d_y(x). First, according to the direction and the magnitude of the edge, the parameter calculator 202 calculates a two-dimensional Gaussian function that is expressed with an anisotropic Gaussian distribution. It is possible to define the two-dimensional Gaussian function as shown in Expressions (3) below, by using a structure tensor H(x) of the derivative value.

$\begin{array}{cc}k\left(x,s\right)=\exp \left(-\frac{1}{h^2}s^TH\left(x\right)s\right)H\left(x\right)=\left[\begin{array}{cc}{d_x\left(x\right)}^2& d_x\left(x\right)d_y\left(x\right)\\ d_x\left(x\right)d_y\left(x\right)& {d_y\left(x\right)}^2\end{array}\right]& \left(3\right)\end{array}$

In Expressions (3), s (where sεN) denotes the position of a point positioned in a predetermined area (hereinafter, a “local neighborhood”) that is centered around the position x. The character “h” (where h>0) denotes a standard deviation of the anisotropic Gaussian distribution. The two-dimensional Gaussian function is strongly affected by the component of the edge in the normal-line direction. As shown in FIG. 4, the two-dimensional Gaussian function has an ellipsoidal shape in which the direction exhibiting a small change in the pixel value is substantially parallel to the long axis, whereas the direction exhibiting a large change in the pixel value is substantially parallel to the short axis. The sharper the edge is, the shorter is the minor axis of the ellipsoidal shape of the two-dimensional Gaussian function (i.e., the ellipse is flattened in the tangent-line direction of the edge). It is possible to calculate the image feature parameters from the structure tensor H(x) by using Expressions (4) below.

$\begin{array}{cc}{\lambda }_{\pm }=\frac{C_{\mathrm{xx}}+C_{\mathrm{yy}}}{2}\pm \sqrt{\frac{{\left(C_{\mathrm{xx}}+C_{\mathrm{yy}}\right)}^2}{4}+C_{\mathrm{xy}}^2}\theta =\{\begin{array}{cc}\frac{\pi }{4}& \mathrm{if}\left(C_{\mathrm{xx}}=C_{\mathrm{yy}}\right)\bigcap \left(C_{\mathrm{xy}}>0\right)\\ -\frac{\pi }{4}& \mathrm{if}\left(C_{\mathrm{xx}}=C_{\mathrm{yy}}\right)\bigcap \left(C_{\mathrm{xy}}<0\right)\\ 0& \mathrm{if}C_{\mathrm{xx}}=C_{\mathrm{yy}}=C_{\mathrm{xy}}=0\\ \frac{1}{2}{\tan }^{-1}\left(\frac{2C_{\mathrm{cy}}}{C_{\mathrm{xx}}-C_{\mathrm{yy}}}\right)& \mathrm{otherwise}\end{array}& \left(4\right)\end{array}$

In this situation, H(x) can be expressed as shown in Expression (5) below.

$\begin{array}{cc}H\left(x\right)=\left[\begin{array}{cc}{d_x\left(x\right)}^2& d_x\left(x\right)d_y\left(x\right)\\ d_x\left(x\right)d_y\left(x\right)& {d_y\left(x\right)}^2\end{array}\right]=\left[\begin{array}{cc}{C}_{\mathrm{xx}}& C_{\mathrm{xy}}\\ C_{\mathrm{xy}}& C_{\mathrm{yy}}\end{array}\right]& \left(5\right)\end{array}$

As shown in FIG. 5, a rotation angle θ is an angle formed by the x axis of the image and the long-axis direction of the two-dimensional Gaussian function. “λ₊” denotes the length of the major axis of the ellipse representing the two-dimensional Gaussian function, whereas “λ₋” denotes the length of the minor axis thereof. It should be noted that “λ₊” and “λ₋” are fixed values for the structure tensor. The long axis of the two-dimensional Gaussian function substantially coincides with the tangent-line direction of the edge. Also, the short axis of the two-dimensional Gaussian function substantially coincides with the normal-line direction of the edge.

It should be noted that by determining any one of a tangent-line direction and a normal-line direction of the edge, the other is determined in a one-to-one correspondence. Therefore, it is possible to use an angle θ′ between an x-axis of the image and a short-axis direction of a two-dimensional Gaussian function instead of the rotation angle θ.

If no further arrangement were made, it would not be possible to calculate the image feature parameters in a stable manner due to the noises contained in the image. Thus, it is acceptable to use a structure tensor obtained by performing a convolution operation with respect to a point positioned within the local neighborhood N that is centered around the position x, as shown in Expression (6) below.

$\begin{array}{cc}H\left(x\right)=\frac{1}{\mathrm{Num}\left(N\right)}\sum _{s\in N}\left[\begin{array}{cc}{d_x\left(x+s\right)}^2& d_x\left(x+s\right)d_y\left(x+s\right)\\ d_x\left(x+s\right)d_y\left(x+s\right)& {d_y\left(x+s\right)}^2\end{array}\right]& \left(6\right)\end{array}$

The local neighborhood N may have an arbitrary shape. For example, it is acceptable to use a rectangular area corresponding to a 5-by-5 tap (i.e., five pixels by five pixels; hereinafter, the notation of “pixels” may be omitted) that is centered around the position x, as the local neighborhood N.

At step S303, according to the rotation angle θ, the converting unit 203 performs a coordinate converting process to convert the coordinates of the position s (hereinafter, the “local coordinates s”) within the local neighborhood N with respect to the position x to local coordinates u expressed with rotated coordinates. It is possible to express the coordinate converting process to convert the x-y coordinates of the image into u-v local coordinates of the two-dimensional Gaussian function by using Expressions (7) below.

$\begin{array}{cc}u=R^{-1}\left(\theta \right)sR^{-1}\left(\theta \right)=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]& \left(7\right)\end{array}$

In this situation, u=(u, v)^T denotes the local coordinates of the two-dimensional Gaussian function.

At step S304, the fitting unit 204 calculates parameters for a curved plane used in a fitting process, by performing a curved-plane fitting process that employs a least-squares method. In other words, the fitting unit 204 performs a fitting process to fit a local gradient around the target pixel into a curved-plane model by the method of least-squares for each of pixels of the input image, and an outputting process to obtain a reconstructed pixel value of the target pixel by using the curved-plane model.

Next, a curved-plane fitting process that employs an n'th-order polynomial curved-plane model, which is disclosed in the reference 1, will be explained. For example, the n'th-order polynomial curved-plane model can be expressed by using Expression (8) shown below.

a₀+a₁s+a₂t+a₃s²+a₄st+a₅t²+ . . .   (8)

In Expression (8), s=(s, t)^T denotes the local coordinates. The curved-plane fitting process using the curved-plane model is performed with respect to the local neighborhood N, while using the position x as the center. An output pixel value obtained after the curved-plane fitting process has been performed can be expressed as I(x)=a₀.

This curved-plane fitting process that is disclosed in the reference 1 has problems as descried below.

For example, let us discuss an edge that extends in a diagonal direction. To perform a fitting process on the edge by using a quadratic (i.e., second-order) curved plane, it is necessary to have all the parameters a₀, a₁, a₂, a₃, a₄, and a₅. In contrast, to perform a fitting process on an edge extending in a longitudinal direction, only the parameters a₀, a₁, and a₃ are necessary, because the t component of the local coordinates is not necessary for the edge extending in the longitudinal direction. Similarly, to perform a fitting process on an edge extending in a transversal direction, only the parameters a₀, a₂, and a₅ are necessary.

As explained here, the fitting process disclosed in the specification of United States Patent Application Publication No. 2007/0047838 (cf. the pamphlet of International Publication No. 2007/027893) is dependent on the direction of each edge. When the number of parameters that are required in the fitting process increases, there is a possibility that the stability of the fitting process may be lowered. In an attempt to cover all the directions, more parameters than necessary will be used, and there is a possibility that over-fitting may occur. The level of stability of the fitting process that employs the least-squares method is determined by the number of parameters to be estimated and the number of sample points that can be fitted. When more parameters than necessary with respect to the number of sample points are used, the fitting process becomes unstable.

The long-axis direction of the two-dimensional Gaussian function corresponds to the tangent-line direction of the edge. The short-axis direction of the two-dimensional Gaussian function corresponds to the normal-line direction of the edge. The change in the pixel value is large in the normal-line direction of the edge and is small in the tangent-line direction of the edge. When this characteristic is applied to the fitting process using a curved plane, the number of parameters for the curved plane becomes large, because the change in the pixel value is large in the normal-line direction of the edge. As a result, the degree of freedom becomes higher, and it is therefore possible to perform the fitting process with a higher level of precision. In contrast, in the tangent-line direction of the edge, the number of parameters for the curved plane does not have to be large, because the change in the pixel value is small. Accordingly, the curved-plane model can be expressed by using Expression (9) shown below.

a₀+a₁v+a₂v²+ . . .   (9)

In this situation, the local coordinates u=(u, v)^T of the two-dimensional Gaussian function are used. The character “u” denotes the long axis of the two-dimensional Gaussian function, whereas the character “v” denotes the short axis of the two-dimensional Gaussian function. In the present example, the degree of freedom of the curved plane is specified only in the short-axis direction of the two-dimensional Gaussian function, while a 0th-order polynomial is used for the long-axis direction. The sharpness of the edge is adjusted in the short-axis direction, while denoising performance (i.e., to eliminate the noises) is improved in the long-axis direction. By using the curved-plane model based on the rotation angle θ, it is possible to realize a fitting process that is not dependent on the direction of the edge. In addition, it is possible to enhance the sharpness of the edge by using fewer parameters than in the example disclosed in the reference 1. Thus, it is possible to improve the level of stability of the fitting process even if the same sample points are used.

For example, according to the technique disclosed in the reference 1, it is necessary to use six parameters to express an edge with a level of precision of a quadratic curved plane. In contrast, according to the first embodiment, only three parameters are necessary. As another example, according to the technique disclosed in the reference 1, it is necessary to use fifteen parameters to express an edge with a level of precision of a quartic (i.e., fourth-order) curved plane. In contrast, according to the first embodiment, only five parameters are necessary. As explained here, with regard to the number of required parameters, the higher the degree of freedom of the curved plane is, the larger is the difference between the conventional technique and the technique according to the first embodiment of the present invention.

The curved-plane fitting process is performed by using the least-squares method. It is possible to define the curved-plane model as shown in Expression (10) below.

$\begin{array}{cc}f\left(u\right)=\left(1vv^2v^3v^4\right)\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ a_4\end{array}\right]=p\left(u\right)a& \left(10\right)\end{array}$

In Expression (10), a quartic polynomial curved-plane model is shown. However, it is acceptable to use a polynomial curved model of any order selected out of 0th to 5th and 7th and higher. Alternatively, instead of using a polynomial, it is acceptable to use a curved-plane model based on a sine wave, as shown in Expression (11) below.

$\begin{array}{cc}f\left(u\right)=\left(1\sin w_0v\sin 2w_0v\sin 3w_0v\sin 4w_0v\right)\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\\ a_4\end{array}\right]=p\left(u\right)a& \left(11\right)\end{array}$

w₀: REFERENCE FREQUENCY

The pixel value of the pixel in the position s within the local neighborhood N centered around the position x will be expressed as I(x+s) (where sεN). The coordinate conversion applied to the local coordinates s according to the rotation angle θ can be realized by performing a coordinate converting step using a rotating matrix while using Expressions (12) shown below.

$\begin{array}{cc}u=R^{-1}\left(\theta \right)sR^{-1}\left(\theta \right)=\left[\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{array}\right]& \left(12\right)\end{array}$

By using Expressions (12), the pixel value I(x+s) and the curved plane f(R⁻¹(θ)s) are brought into correspondence with each other. The least-squares method is a method for obtaining a parameter that makes the squared error in this correspondence the smallest and can be defined as shown in Expressions (13) below.

$\begin{array}{cc}\begin{array}{c}\hat{a}\left(x\right)=\underset{a}{\min }E\left(x,a\right)\\ E\left(x,a\right)=\sum _{s\in N}k\left(x,s\right){I\left(x+s\right)-f\left(R^{-1}\left(\theta \right)s\right)}^2\\ =\sum _{s\in N}k\left(x,s\right){I\left(x+s\right)-p\left(R^{-1}\left(\theta \right)s\right)a}^2\end{array}& \left(13\right)\end{array}$

• • â(x) IS WRITTEN AS â(x) IN THE SPECIFICATION

In Expressions (13), a â(x) is a parameter used in the fitting process that employs the least-squares method. k(x, s) denotes a weight applied to the point s. In the present example, a two-dimensional Gaussian function is used. The local neighborhood N may have any arbitrary shape. For example, it is acceptable to use a rectangular area corresponding to a 5-by-5 tap that is centered around the position x, as the local neighborhood N. In the present example, to simplify the explanation, Expressions (13) will be expressed in a matrix format as shown in Expressions (14) below.

$\begin{array}{cc}E\left(x,a\right)={\left(Y-Pa\right)}^TW\left(Y-Pa\right)Y=\left[\begin{array}{c}I\left(x+s_0\right)\\ ⋮\\ I\left(x+s_n\right)\end{array}\right]W=\left[\begin{array}{ccc}k\left(x,s_0\right)& 0& 0\\ 0& \ddots & 0\\ 0& 0& k\left(x,s_n\right)\end{array}\right]P=\left[\begin{array}{c}p\left(R^{-1}\left(\theta \right)s_0\right)\\ ⋮\\ p\left(R^{-1}\left(\theta \right)s_n\right)\end{array}\right]& \left(14\right)\end{array}$

In this situation, points within the local neighborhood N are expressed as N={s₀, . . . , s_n}. When the matrix format is used, it is possible to find a solution to the least-squares method in a one-to-one correspondence by performing a matrix calculation as shown in Expression (15) below.

â(x)=(P^TWP)⁻¹P^TWY  (15)

Expression (15) is called a normal equation. When a linear least-squares method is used, the normal equation provides the optimal solution. It is possible to calculate values for an inverse matrix by using, for example, a Lower Upper (LU) decomposition or a singular value decomposition. Let us assume that a â(x) has been calculated as shown in Expression (16) below.

â(x)=(â₀(x),â₁(x),â₂(x),â₃(x),â₄(x))^T  (16)

The output pixel obtained after the fitting process has been performed can be calculated by using Expression (17) shown below.

â₀(x)  (17)

Next, operations performed in the curved-plane fitting step at step S304 will be explained, with reference to FIG. 6. At step S601, the fitting unit 204 calculates a matrix P according to Expressions (14). At step S602, the fitting unit 204 calculates a matrix W according to Expressions (14) by using the image feature parameters shown in Expressions (4). At step S603, the fitting unit 204 calculates (P^TWP)⁻¹ by using an LU decomposition, a singular value decomposition, or the like. At step S604, the fitting unit 204 calculates a â(x) according to Expression (15).

Next, a result of a fitting process performed by using the conventional kernel regression method and a result of the fitting process according to the first embodiment of the present invention will be explained with reference to FIGS. 7A, 7B, and 8A to 8D.

FIGS. 8A and 8B are drawings illustrating results obtained by performing a curved-plane fitting process on the image shown in FIG. 7B by using the least-squares method and performing a first-order derivative operation in the x direction and in the y direction, respectively. By referring to FIGS. 8A and 8B, it is understood that the curve parameters are distributed to the axes, while the model parameters are dependent on the rotation angle of the two-dimensional Gaussian function. As a result, to reproduce the edges in all the directions at the same level, it is necessary to have all the parameters on the quadratic curved plane.

FIGS. 8C and 8D are drawings illustrating results obtained by performing a curved-plane fitting process on the image shown in FIG. 7B by using the least-squares method and performing a first-order derivative operation in the u direction (i.e., the long-axis direction of the two-dimensional Gaussian function) and in the v direction (i.e., the short-axis direction of the two-dimensional Gaussian function), respectively. By referring to FIGS. 8C and 8D, it is understood that the axes in the curved-plane model coincide with the main component axes of the edges. As a result, the components of the edges are concentrated on the v axis (i.e., the short axis of the two-dimensional Gaussian function). Accordingly, to reproduce the edges in all the directions, using only the v axis is sufficient.

Next, differences between the conventional method and the method according to the first embodiment of the present invention will be explained, with reference to FIG. 9. According to the conventional method, the fitting process based on kernel regression is performed on the x-y axes. In contrast, according to the first embodiment of the present invention, the fitting process based on kernel regression is performed on the u-v axes. Because the u-v axes are the main component axes of the edges, it is possible to distribute the degree of freedom in the maximum-amplitude direction (i.e., the main components) of the information.

As explained above, according to the first embodiment, it is possible to effectively eliminate noises from the image while maintaining the shapes and the edges in the image.

Next, a second embodiment of the present invention will be explained. Some of the constituent elements of an image processing unit 102b according to the second embodiment shown in FIG. 10 that are the same as those of the image processing unit 102a shown in FIG. 2 will be referred to by using the same reference characters. The second embodiment is different from the first embodiment in that the image processing unit 102b includes a selecting unit 1001.

According to the first embodiment, the filtering process that is excellent in terms of reproducibility of the edges is performed by rotating the curved-plane model in the tangent-line directions of the edges. The reason for this can be explained as follows: In the edge portions, the information is concentrated in the short-axis direction of the ellipse of the two-dimensional Gaussian function (i.e., the normal-line direction of the edge) as explained above. As a result, the curved-plane in which the degree of freedom is given to the short-axis component is fitted well.

It should be considered, however, that an image contains not only edge areas, but also flat areas and corner areas (e.g., corners and pointed extremities). In such areas, the information is not necessarily concentrated on the short-axis component. According to the second embodiment, areas are classified according to image features, so that an appropriate one of curved-plane models can be assigned according to the classification. According to the second embodiment, the image features are obtained by using an x-direction derivative value and a y-direction derivative value. More specifically, the image features are obtained by using a structure tensor, which is also explained in the description of the first embodiment.

As shown in FIG. 11, the operations performed by the image processing unit 102b shown in FIG. 10 are different from the operations shown in FIG. 3 in, at least, that a curved-plane model selecting step at step S1103 is performed before the coordinate converting step and the curved-plane fitting step. The other steps will be explained while a focus is placed on the differences from FIG. 3.

At step S1102, the tangent-line direction of the edge as well as the parameters corresponding to the long axis and the short axis of the ellipse that express local features of the image are calculated by using a structure tensor. Harris et al. have presented a feature classification of images based on a structure tensor (cf. C. Harris and M. Stephens (1988), “A Combined Corner and Edge Detector”, Proceedings of the 4th Alvey Vision Conference, pp. 147-151).

According to the feature classification, when fixed values of the structure tensor are expressed as λ₊ and λ₋, each of the portions in an image can be classified as an edge area (“Edge”), a flat area (“Flat”), or a corner area (“Corner”) as shown in FIG. 12, according to the fixed values λ₊ and λ₋. When one of the fixed values λ₊ and λ₋ is large and the other is small, the ellipse of the two-dimensional Gaussian function has a flattened shape. In that situation, the corresponding image portion is classified as an edge area. When both of the fixed values are large, the ellipse of the two-dimensional Gaussian function has an isotropic circular shape and is small. In that situation, the corresponding image portion is classified as a corner area (including a corner or a pointed extremity). When both of the fixed values are small, the ellipse of the two-dimensional Gaussian function has an isotropic circular shape and is large. In that situation, the corresponding image portion is classified as a flat area.

The curved-plane model according to the first embodiment is a curved-plane model that uses only the short-axis component of the ellipse of the two-dimensional Gaussian function. This curved-plane model is based on the notion that the information is concentrated in the normal-line directions of the edges in the edge areas. Thus, the model is suitable for the edge areas among the areas explained above. For example, according to the first embodiment, an edge area curved-plane model as shown in Expression (18) below may be used.

f_Edge(u)=a₀+a₁v+a₂v²+a₃v³+a₄v⁴  (18)

In contrast, in corner areas and flat areas, the ellipse of the two-dimensional Gaussian function is isotropic, and the information is not necessarily concentrated on the short axis (i.e., the v axis). Thus, it is more appropriate to use both the u axis and the v axis in the corner areas and the flat areas. Of these areas, with regard to the corner areas, because the change in the pixel value is large, it is more appropriate to give a higher degree of freedom to the curved-plane model. For example, a corner area curved-plane model as shown in Expression (19) below may be used.

f_corner(u)=a₀+a₁u+a₂u²+a₃u³+a₄u⁴+a₅v+a₆v²+a₇v³+a₈v⁸  (19)

In contrast, because it is desirable to keep the noises as little as possible in the flat areas, a flat area curved-plane model based on a low-order polynomial as shown in Expression (20) below may be used.

f_Flat(u)=a₀+a₁u+a₂u²+a₃v+a₄v²  (20)

At step S1103, the selecting unit 1001 classifies a target pixel as one of a plurality of predetermined classes (i.e., areas). The selecting unit 1001 further selects one of the curved-plane models that respectively correspond, in a one-to-one correspondence, to the plurality of classes used in the classification, as the curved-plane model for the target pixel. More specifically, the selecting unit 1001 classifies the target image based on the fixed values λ₊ and λ₋ of the structure tensor. The classification process is performed according to, for example, the classification explained with reference to FIG. 12. The selecting unit 1001 selects the curved-plane model as explained above according to the result of the classification process.

At step S1105, the fitting unit 204 performs the fitting process that employs the least-squares method, by using the curved-plane model that has been selected at step S1103.

Next, a third embodiment of the present invention will be explained. Some of the constituent elements of an image processing unit 102c according to the third embodiment shown in FIG. 13 that are the same as those of the image processing unit 102a shown in FIG. 2 will be referred to by using the same reference characters. The third embodiment is different from the first embodiment in that the image processing unit 102c includes a filter selecting unit 1301 and a convolution operating unit 1302, instead of the converting unit 203 and the fitting unit 204.

According to the first embodiment, it is necessary to solve the normal equation for each of the pixels that are the targets of the filtering process. More specifically, it is necessary to calculate the values for the inverse matrix by using an LU decomposition or a singular value decomposition. According to the third embodiment, for the purpose of reducing the processing cost, results obtained by solving the normal equation in advance are stored in a Look-Up Table (LUT). With this arrangement, it is possible to enhance viability in a circuit or the like.

The operations shown in FIG. 14 that are performed by the image processing unit 102c shown in FIG. 13 are different from the operations explained with reference to FIG. 3 in, at least, that the procedure includes a filter selecting step at step S1403 and a convolution operating step at step S1404, instead of the coordinate converting step at step S303 and the curved-plane fitting step at step S304. The other steps are the same as those shown in FIG. 3.

At step S1403, the filter selecting unit 1301 selects an appropriate filter out of the LUT showing the results obtained by solving the normal equation, based on the image feature parameters that have been calculated at step S1402.

The normal equation is obtained by using Expression (15). As shown in Expressions (14), Y denotes a pixel value and varies depending on the input image. The value (P^rWP)⁻¹P^rW is dependent only on the image feature parameters λ₊, λ₋, and θ and is not dependent on the image. The two-dimensional Gaussian function can be expressed by using self-correlation coefficients for a derivative of the input image, as shown in Expressions (21).

$\begin{array}{cc}k\left(x,s\right)=\exp \left(-\frac{1}{h^2}s^TH\left(x\right)s\right)H\left(x\right)=\left[\begin{array}{cc}{d_x\left(x\right)}^2& d_x\left(x\right)d_y\left(x\right)\\ d_x\left(x\right)d_y\left(x\right)& {d_y\left(x\right)}^2\end{array}\right]& \left(21\right)\end{array}$

By rewriting Expressions (21) with the image feature parameters, Expression (22) can be obtained.

$\begin{array}{cc}\begin{array}{c}H\left(x\right)=\left[\begin{array}{cc}{\lambda }_{+}{\cos }^2\theta +{\lambda }_{-}{\sin }^2\theta & \left({\lambda }_{+}-{\lambda }_{-}\right)\sin \theta \cos \theta \\ \left({\lambda }_{+}-{\lambda }_{-}\right)\sin \theta \cos \theta & {\lambda }_{+}{\sin }^2\theta +{\lambda }_{-}{\cos }^2\theta \end{array}\right]\\ =H\left({\lambda }_{+},{\lambda }_{-},\theta \right)\end{array}& \left(22\right)\end{array}$

It is also possible to rewrite the two-dimensional Gaussian function as shown in Expression (23) below.

$\begin{array}{cc}k\left({\lambda }_{+},{\lambda }_{-},\theta ,s\right)=\exp \left(-\frac{1}{h^2}s^TH\left({\lambda }_{+},{\lambda }_{-},\theta \right)s\right)& \left(23\right)\end{array}$

It is possible to express the matrix W as shown in Expression (24) below.

$\begin{array}{cc}\begin{array}{c}W\left({\lambda }_{+},{\lambda }_{-},\theta \right)=\left[\begin{array}{ccc}k\left(x,s_0\right)& 0& 0\\ 0& \ddots & 0\\ 0& 0& k\left(x,s_n\right)\end{array}\right]\\ =\left[\begin{array}{ccc}k\left({\lambda }_{+},{\lambda }_{-},\theta ,s_0\right)& 0& 0\\ 0& \ddots & 0\\ 0& 0& k\left({\lambda }_{+},{\lambda }_{-},\theta ,s_n\right)\end{array}\right]\end{array}& \left(24\right)\end{array}$

It is understood that the matrix W is dependent only on the image feature parameters. Similarly, it is possible to express the matrix P as shown in Expression (25) below.

$\begin{array}{cc}P\left(\theta \right)=\left[\begin{array}{c}p\left(R^{-1}\left(\theta \right)s_0\right)\\ ⋮\\ p\left(R^{-1}\left(\theta \right)s_n\right)\end{array}\right]& \left(25\right)\end{array}$

It is understood that the matrix P is dependent only on the image feature parameters. Accordingly, it is understood that Expression (26) is also dependent only on the image feature parameters.

X(λ₊,λ₋,θ)=(P(θ)^TW(λ₊,λ₋,θ)P(θ))⁻¹P(θ)^TW(λ₊,λ₋,θ)  (26)

By calculating X(λ₊, λ₋, θ)₁ (where 1=0, . . . , L), in advance, with respect to an arbitrary set of image feature parameters that have been discretized such as (λ₊, λ₋, θ)₁ (where 1=0, . . . , L), it is possible to obtain the solution by using Expression (27) without performing any additional calculation.

â(x)=X(λ₊,λ₋,θ)₁Y  (27)

In other words, at the filter selecting step at step S1403, the values X(λ₊, λ₋, θ)₁ (where 1=0, . . . , L) are calculated in advance with respect to (λ₊, λ₋, θ)₁ (where 1=0, . . . , L) so that the calculated values are registered into the LUT (i.e., a filter bank). In the actual process, a corresponding one of the values X(λ₊, λ₋, θ)₁ is selected out of the LUT, based on the calculated image feature parameters (λ₊, λ₋, θ).

At step S1404, the convolution operating unit 1302 performs a convolution operation with the pixel value vector Y, by using the filter X(λ₊, λ₋, θ)₁ that has been selected at the filter selecting step at step S1403. Further, the convolution operating unit 1302 calculates an output pixel to which a curved plane has been fitted by using the least-squares method. More specifically, the convolution operating unit 1302 performs a matrix calculation shown in Expression (28) below.

$\begin{array}{cc}\begin{array}{c}\hat{a}\left(x\right)={X\left({\lambda }_{+},{\lambda }_{-},\theta \right)}_1Y\\ ={X\left({\lambda }_{+},{\lambda }_{-},\theta \right)}_1\left[\begin{array}{c}I\left(x+s_0\right)\\ ⋮\\ I\left(x+s_n\right)\end{array}\right]\end{array}& \left(28\right)\end{array}$

Next, a fourth embodiment of the present invention will be explained. Some of the constituent elements of an image processing unit 102d according to the fourth embodiment shown in FIG. 15 that are the same as those of the image processing unit 102a shown in FIG. 2 will be referred to by using the same reference characters. The fourth embodiment is different from the first embodiment in that the image processing unit 102d includes a weight calculator 1501.

According to the first and the second embodiments, the filtering process that is suitable for the edges in the input image is performed by using the two-dimensional Gaussian function. However, because it is not possible to express, for example, apex portions of a triangle by using a two-dimensional Gaussian function, the image becomes dull. To cope with this problem, according to the fourth embodiment, the two-dimensional Gaussian function is made robust based on pixel values of the pixels that are in the apex portions of a triangle or the like. As a result, it is possible to solve the problem where, for example, the apex portions become dull.

As shown in FIG. 16, the operations performed by the image processing unit 102d shown in FIG. 15 are different from the operations explained with reference to FIG. 3 in, at least, that the procedure includes a weight calculating step at step S1606 after the curved-plane fitting step at step S1605. The other steps will be explained while a focus is placed on the differences from FIG. 3.

The process performed at step S1605 is basically the same as the process performed according to the first embodiment, except that the fitting unit 204 calculates the matrix W by using Expression (29) shown below.

$\begin{array}{cc}W=\left[\begin{array}{ccc}w\left(s_0\right)k\left(x,s_0\right)& 0& 0\\ 0& \ddots & 0\\ 0& 0& w\left(s_n\right)k\left(x,s_n\right)\end{array}\right]& \left(29\right)\end{array}$

At step S1606, the weight calculator 1501 applies a weight to a pixel value in the input image within the local neighborhood, according to the distance from the target pixel that is positioned at the center of the local coordinates. For example, the weight calculator 1501 calculates the weight in the local neighborhood in such a manner that the larger the difference from the pixel value of the target pixel in the input image is, the smaller is the weight (i.e., the smaller the difference is, the larger is the weight). Two-dimensional Gaussian functions are applicable to step edge portions in an image. However, two-dimensional Gaussian functions may not be applicable to extremity portions and corners in some situations. When a fitting process that employs the least-squares method is performed in such areas to which two-dimensional Gaussian functions are not applicable, the image is severely distorted because of out-of-range pixel values to which the function is not applicable. As a result, a problem arises where, for example, corners become dull. The areas having such out-of-range values are assumed to be areas containing pixels each having a pixel value that is different from the pixel value of the target pixel. It is possible to define the weight based on the difference in the pixel values by using a Gaussian function as shown in Expressions (30) below.

$\begin{array}{cc}w\left(s\right)=k_{\sigma 1}\left(I\left(x+s\right)-I\left(s\right)\right)k_{\sigma 1}\left(I\right)=\exp \left(-\frac{I^TI}{{\sigma }_1^2}\right)& \left(30\right)\end{array}$

In this situation, σ_I (where σ_I>0) denotes a standard deviation of the Gaussian function and serves as a parameter indicating what degree of difference in the pixel values causes the weight to be smaller.

As another example, in a situation where a block distortion occurs (e.g., in an encoded image), the weight calculated by using Expressions (30) also reacts to the block distortion. As a result, it is not possible to solve the problem of the block distortion. To cope with this situation, it is also acceptable to use a deblocking-type weight as shown in Expression (31) below so that the weight based on the difference in the pixel values is not used in block boundary areas.

$\begin{array}{cc}w\left(s\right)=\{\begin{array}{cc}{k}_{\sigma 1}\left(I\left(x+s\right)-I\left(s\right)\right)& \begin{array}{c}\mathrm{WHEN}x\mathrm{AND}s\mathrm{ARE}\\ \begin{array}{c}\mathrm{POINTS}\mathrm{IN}\mathrm{THE}\\ \mathrm{SAME}\mathrm{BLOCK}\end{array}\end{array}\\ 1& \mathrm{OTHERWISE}\end{array}& \left(31\right)\end{array}$

Next, a fifth embodiment of the present invention will be explained. Some of the constituent elements of an image processing unit 102e according to the fifth embodiment shown in FIG. 17 that are the same as those of the image processing unit 102a shown in FIG. 2 will be referred to by using the same reference characters. The fifth embodiment is different from the first embodiment in that the image processing unit 102e includes a difference calculator 1701.

According to the fifth embodiment, a difference between a pixel value on the curved plane that is fitted according to any of the first to the fourth embodiments and a pixel value in the input image is calculated so that a fitting process is performed again on the calculated difference. With this arrangement, it is possible to enhance reproducibility of, for example, a fine texture.

As shown in FIG. 18, the operations performed by the image processing unit 102e shown in FIG. 17 are different from the operations explained with reference to FIG. 3 in, at least, that the procedure includes a difference calculating step at step S1806, a coordinate calculating step at step S1807, and a curved-plane fitting step at step S1808, after the curved-plane fitting step at step S1805. The other steps will be explained while a focus is placed on the differences from FIG. 3.

In this situation, when the input image is assumed to be made up of a framework component, a texture component, and noises, it is possible to provide a model for the pixel I(x) as shown in Expression (32) below.

I(x)=S(x)+T(x)+n(x)  (32)

In Expression (32), S(x) denotes a framework component, whereas T(x) denotes a texture component, while n(x) denotes a noise. Expression (32) corresponds to an addition model. It is also acceptable to use a multiplication model as shown in Expression (33) below.

I(x)=S(x)·T(x)+n(x)  (33)

According to the first to the fourth embodiments, the degree of freedom of the curved plane is specified in the normal-line directions of the edges and is not specified in the tangent-line directions of the edges. This arrangement is equivalent to extracting the framework component while ignoring irregularities in the tangent-line directions of the edges. Thus, it is possible to consider the result of the curved-plane fitting process that employs the least-squares method as the framework component of the image as shown in Expression (34) below.

S(x)=â₀(x)  (34)

At step S1806, the difference calculator 1701 calculates the difference between a pixel value in the input image and a pixel value of the framework component. First, the texture component can be calculated as shown in Expression (35) below, based on Expressions (32) and (34).

T(x)=I(x)−â₀(x)−n(x)  (35)

From this, it is understood that, to obtain the texture component, it is necessary to obtain the difference between the pixel value in the input image and the pixel value of the framework component. A fitting curved plane can be expressed as shown in Expression (36) below.

f(u)=p(u)a  (36)

Accordingly, it is possible to calculate the difference by using Expression (37) below.

$\begin{array}{cc}\begin{array}{c}r\left(x+s\right)=I\left(x+s\right)-f\left(R^{-1}\left(\theta \right)s\right)\\ =I\left(x+s\right)-p\left(R^{-1}\left(\theta \right)s\right)\hat{a}\end{array}& \left(37\right)\end{array}$

In the present example, s denotes a point within the local neighborhood N.

At step S1808, the fitting unit 204 calculates the texture component by performing, with respect to the difference that has been calculated at step S1806, a curved-plane fitting process that employs the least-squares method on the local coordinates obtained as a result of the coordinate converting process performed at step S1807. It is possible to express a texture curved-plane model by using Expression (38) shown below.

$\begin{array}{cc}{f}_r\left(u\right)=\left(1uu^2u^3u^4\right)\left[\begin{array}{c}{b}_0\\ b_1\\ b_2\\ b_3\\ b_4\end{array}\right]=p_r\left(u\right)b& \left(38\right)\end{array}$

According to the first to the fourth embodiments, the model is specified in the normal-line directions of the edges. In contrast, in the present example, the model is specified in the tangent-line directions of the edges. This arrangement is made so that it is possible to extract the texture component in the tangent-line directions of the edges that are orthogonal to the normal-line directions of the edges based on which the framework component of the image is extracted. It is possible to express the least-squares method applied to the difference by using Expressions (39) shown below.

$\begin{array}{cc}\begin{array}{c}\hat{b}\left(x\right)=\underset{b}{\min }E\left(x,b\right)\\ E\left(x,b\right)=\sum _{s\in N}k\left(x+s\right){r\left(x,s\right)-f_r\left(R^{-1}\left(\theta \right)s\right)}^2\\ =\sum _{s\in N}k\left(x+s\right){r\left(x,s\right)-p_r\left(R^{-1}\left(\theta \right)s\right)b}^2\end{array}& \left(39\right)\end{array}$

• • {circumflex over (b)}(x) IS WRITTEN AS b̂(x) IN THE SPECIFICATION

By putting Expressions (39) into a matrix format, Expressions (40) below can be obtained.

$\begin{array}{cc}E\left(x,b\right)={\left(Y_r-P_rb\right)}^TW\left(Y_r-P_rb\right)Y_r=\left[\begin{array}{c}r\left(x+s_0\right)\\ ⋮\\ r\left(x+s_n\right)\end{array}\right]W=\left[\begin{array}{ccc}k\left(x,s_0\right)& 0& 0\\ 0& \ddots & 0\\ 0& 0& k\left(x,s_n\right)\end{array}\right]P_r=\left[\begin{array}{c}{p}_r\left(R^{-1}\left(\theta \right)s_0\right)\\ ⋮\\ p_r\left(R^{-1}\left(\theta \right)s_n\right)\end{array}\right]& \left(40\right)\end{array}$

Accordingly, it is possible to express a normal equation as shown in Expression (41) below.

{circumflex over (b)}(x)=P_r^TWP_r)⁻¹P_r^TWY_r  (41)

It is possible to calculate values for an inverse matrix by using, for example, an LU decomposition or a singular value decomposition. Let us assume that b̂(x) has been calculated as shown in Expression (42) below.

{circumflex over (b)}(x)=({circumflex over (b)}₀(x),{circumflex over (b)}₁(x),{circumflex over (b)}₂(x),{circumflex over (b)}₃(x),{circumflex over (b)}₄(x))  (42)

In this situation, the texture component can be expressed by using Expression (43) below.

{circumflex over (b)}₀(x)  (43)

The image processing unit 102e adds the curved-plane model and the texture curved-plane model together and generates an output image. In other words, it is possible to express the output image that is obtained after the filtering process has been performed, by using Expression (44) below.

â₀(x)+{umlaut over (b)}₀(x)  (44)

The present invention is not limited to the exemplary embodiments described above. It is possible to apply modifications to the present invention without departing from the gist thereof. An image processing computer program used by any of the image processing apparatuses according to the embodiments described above to execute the processing steps may be provided as a computer program product as being recorded on a computer-readable recording medium such as a Compact Disk Read-Only Memory (CD-ROM), a Flexible Disk (FD), a Compact Disk Recordable (CD-R), or a Digital Versatile Disk (DVD), in a file that is in an installable format or in an executable format.

Further, another arrangement is acceptable in which the image processing computer program according to any of the embodiments described above is provided as being incorporated in a ROM or the like in advance.

The image processing computer program executed by any of the image processing apparatuses according to the embodiments described above has a module configuration that includes the image processing unit 102 described above. As the actual hardware configuration, the image processing unit 102 is loaded into a main storage device (not shown) when a Central Processing Unit (CPU) (i.e., a processor; not shown) reads and executes the image processing computer program from a storage medium (not shown), so that the image processing unit 102 is generated in the main storage device.

Additional advantages and modifications will readily occur to those skilled in the art. Therefore, the invention in its broader aspects is not limited to the specific details and representative embodiments shown and described herein. Accordingly, various modifications may be made without departing from the spirit or scope of the general inventive concept as defined by the appended claims and their equivalents.

Claims

1. An image processing apparatus comprising:

a calculating unit that calculates at least one of a tangent-line direction and a normal-line direction of an edge around a target pixel in an input image by using a direction and a magnitude of a gradient of pixel values of pixels that neighbor the target pixel in the input image;

a converting unit that converts coordinates of pixels around the target pixel into a rotated coordinate system by rotating the axes according to an angle between a horizontal direction of the input image and the tangent-line direction of the edge or the normal-line direction of the edge; and

a fitting unit that performs a fitting process to fit a local gradient around the target pixel, on the rotated coordinate system, into a curved-plane model by the method of least-squares, and an outputting process to obtain a reconstructed pixel value of the target pixel by using the curved-plane model.

2. The apparatus according to claim 1, further comprising:

a selecting unit that classifies the target pixel as one of a plurality of predetermined classes by using the direction and the magnitude of the gradient, and selects one of curved-plane models that respectively correspond, in a one-to-one correspondence, to the classes as a curved-plane model for the target pixel, wherein

the fitting unit performs the fitting process by using the curved-plane model selected by the selecting unit.

3. The apparatus according to claim 1, wherein the fitting unit configures the curved-plane model so that a degree of freedom in a normal-line direction of the edge is higher than a degree of freedom in a tangent-line direction of the edge.

4. The apparatus according to claim 1, further comprising a weight calculator that applies a weight to a pixel value of a pixel in the input image according a distance from the target pixel to the pixel.

5. The apparatus according to claim 4, wherein the weight calculator calculates the weight in such a manner that the larger a difference between the pixel value of the pixel in the input image and the pixel value of the target pixel is, the smaller is the weight applied to the pixel value of the pixel in the input image.

6. The apparatus according to claim 1, further comprising:

a storage unit that stores a normal equation based on the least-squares method for each of mutually-different tangent-line directions of the edge, wherein

the fitting unit performs the fitting process using the curved-plane model by selecting a normal equation corresponding to a tangent-line direction of the edge from the normal equations stored in the storage unit and performing a convolution operation on the selected normal equation and the pixel value of the target pixel at the rotated coordinates.

7. The apparatus according to claim 1, further comprising:

a difference calculator that calculates a difference between the pixel value of the pixel in the input image and a pixel value in the curved-plane model, wherein

the fitting unit performs the fitting process that employs the least-squares method by using a texture curved-plane model that is formed according to the pixel value of the target pixel, based on the calculated difference.

8. An image processing method comprising:

calculating at least one of a tangent-line direction and a normal-line direction of an edge around a target pixel in an input image by using a direction and a magnitude of a gradient of pixel values of pixels that neighbor the target pixel in the input image;

converting coordinates of pixels around the target pixel into a rotated coordinate system by rotating the axes according to an angle between a horizontal direction of the input image and the tangent-line direction of the edge or the normal-line direction of the edge; and

performing a fitting process to fit a local gradient around the target pixel, on the rotated coordinate system, into a curved-plane model by the method of least squares, and an outputting process to obtain a reconstructed pixel value of the target pixel by using the curved-plane model.

9. A computer program product having a computer readable medium including programmed instructions for performing image processing, wherein the instructions, when executed by a computer, cause the computer to perform:

calculating at least one of a tangent-line direction and a normal-line direction of an edge around a target pixel in an input image by using a direction and a magnitude of a gradient of pixel values of pixels that neighbor the target pixel in the input image;

converting coordinates of pixels around the target pixel into a rotated coordinate system by rotating the axes according to an angle between a horizontal direction of the input image and the tangent-line direction of the edge or the normal-line direction of the edge; and

performing a fitting process to fit a local gradient around the target pixel, on the rotated coordinate system, into a curved-plane model by the method of least squares, and an outputting process to obtain a reconstructed pixel value of the target pixel by using the curved-plane model.