WO2022190249A1 - Signal processing device, signal processing method, and signal processing program - Google Patents

Abstract

This signal processing device is provided with a control unit for generating, from an input signal, an output signal that has undergone edge-preserving smoothing by: applying, to the input signal, a sparse regularization function which is a nonconvex function and which changes the smoothness of gradation by adjusting a first control parameter included in weights of the nonconvex function and changes the steepness of an edge by adjusting a second control parameter included in the weights of the nonconvex function; and calculating a solution to an optimization problem using global optimization.

Description

SIGNAL PROCESSING DEVICE, SIGNAL PROCESSING METHOD AND SIGNAL PROCESSING PROGRAM

The present invention relates to a signal processing device, a signal processing method, and a signal processing program.

Edge-preserving smoothing is one of the important processes in image processing. Edge-preserving smoothing captures edges and contours to reveal image structure and smooth out unimportant details. Edge-preserving smoothing not only improves image quality, but is also used as preprocessing for various image processing such as contour extraction, segmentation, and style conversion. Therefore, various edge-preserving smoothing techniques have been proposed so far, as shown in FIG.

As shown in FIG. 9, the edge-preserving smoothing technique is classified into a local filter technique and a global optimization technique. Global optimization methods use the statistical properties of the entire image to formulate an optimization problem, find a solution, and calculate the output to perform edge-preserving smoothing. Although the global optimization method requires a lot of execution time, it has high edge preserving performance, less artifacts, and has better smoothing performance than the local filter method.

The global optimization approach uses several control parameters. Many techniques of global optimization adjust the smoothness of an image, for example, by changing the value of a control parameter. For example, the control parameter λ of Total Variation (hereinafter also referred to as “TV”) balances the fidelity term and the TV term to adjust the smoothness of the gradation (see Non-Patent Documents 1 and 2, for example). On the other hand, in the _L0 gradient projection, by changing the control parameter α, the total number of pixels in which edges exist is changed to adjust sharp edges (see, for example, Non-Patent Document 3).

However, there is no global optimization method that adjusts both gradation smoothness and edge steepness. Therefore, if the smoothness of the gradation is strengthened using a global optimization method that adjusts the smoothness of the gradation, there is a problem that the edges become dull and the output becomes blurred. On the other hand, if the sharpness of the edge is strengthened by using the global optimization method for adjusting the steepness of the edge, there is a problem that the gradation deteriorates into a false contour.

In view of the above circumstances, the present invention aims to provide a technology that can independently adjust both the smoothness of gradation and the steepness of edges in edge-preserving smoothing by global optimization.

According to one aspect of the present invention, the smoothness of the gradation is changed by adjusting a first control parameter included in the weight of the non-convex function, which is a sparse regularization function of the non-convex function for the input signal. , by adjusting a second control parameter contained in the weight of said non-convex function, applying a sparse regularization function of the non-convex function that varies the steepness of the edges to solve the optimization problem by global optimization. A signal processing apparatus comprising a control unit for generating an edge-preserving smoothed output signal from the input signal by calculating

One aspect of the present invention is a signal processing method performed by a signal processing device, which is a sparse regularization function that is a non-convex function for an input signal, and a first control parameter included in the weight of the non-convex function is applying a non-convex sparse regularization function that varies the smoothness of the gradation by adjusting and the steepness of the edges by adjusting a second control parameter contained in the weights of said non-convex function. and calculating a solution of an optimization problem by global optimization to generate an edge-preserving smoothed output signal from the input signal.

One aspect of the present invention is a signal processing program for causing a computer to execute the above signal processing device.

The present invention makes it possible to independently adjust both gradation smoothness and edge steepness in edge-preserving smoothing by global optimization.

FIG. 4 is a diagram comparing contour lines of the cusp norm according to the present embodiment with contour lines of the L ₁ norm and the L _p pseudo-norm, which are examples of other norms. 2 is a block diagram showing the internal configuration of the signal processing device according to the embodiment; FIG. 4 is a block diagram showing the internal configuration of a first variable updating unit according to the embodiment; FIG. It is a block diagram which shows the internal structure of the 2nd variable update part by this embodiment. It is a block diagram which shows the internal structure of the 3rd variable update part by this embodiment. 4 is a flow chart showing the flow of processing by the signal processing device according to the embodiment; It is a figure (part 1) which shows the implementation result using the signal processing apparatus of this embodiment. It is a figure (part 2) which shows the implementation result using the signal processing apparatus of this embodiment. FIG. 10 illustrates types of hedge-preserving smoothing methods;

Hereinafter, embodiments of the present invention will be described with reference to the drawings. As described above, the present invention aims to independently adjust both gradation smoothness and edge steepness in edge-preserving smoothing by global optimization. To this end, we propose to use the cusp norm, which will be described later, to allow us to incorporate two parameters, one for adjusting the smoothness of the gradation and the other for adjusting the sharpness of the edge. is doing. In addition, since the cusp norm is a non-convex function, it generally takes a lot of computation time to solve the optimization problem formulated in edge-preserving smoothing by global optimization. is known, and a technique for solving this problem of computation time is proposed below. First, the optimization problem formulated in the method of edge-preserving smoothing by global optimization used in this embodiment and its solution will be described. An input color image signal having three channels of RGB is defined by the following equation (1) as an object of hedge preserving smoothing.

In equation (1), N=N _v ×N _h , where N _v is the number of pixels in the vertical direction and N _h is the number of pixels in the horizontal direction. Since there are three channels of red, green, and blue, the formula (1) is "3N". In equation (1), the vector f is a column vector obtained by vertically scanning the input color image signal for each channel and accumulating them, and is hereinafter referred to as the input image signal vector f. An output image signal vector x obtained by performing hedge preserving smoothing of the present embodiment on an input image signal vector f is defined by the following equation (2).

The following equation (3) is an equation representing the edge-preserving smoothing optimization problem of this embodiment expressed using the input image signal vector f and the output image signal vector x.

The function of the second term in Equation (3) is the pointed star gradient (hereinafter also referred to as "PSG" (Pointed Star Gradient)). The name "pointed star" is derived from the pointed star norm, which will be described later. The subscript vector w of the function is a vector containing N weight parameters defined by the following equation (4). The vector w is hereinafter referred to as a weight vector w.

Here, the elements included in the output image signal vector x, that is, the pixel values can be expressed as the following equation (5).

In equation (5), the subscript “i” is the pixel position in the vertical direction, the subscript “j” is the pixel position in the horizontal direction, and the subscript “c” is the color, That is, red, blue, or green.

The function of the second term of formula (3) is defined as the following formula (6).

In Expression (6), the function Grad(·) is a function that calculates the gradient of each pixel of the image signal vector given as an argument. The vec operator is an operator that vertically arranges each column of a matrix given as an argument. When the output image signal vector x of the formula (5) is given as an argument to the function Grad(·), the output of the function Grad(·) at the pixel positions “i” and “j” is given by the following formula (7 ).

However, if i+1>N _v in Equation (7), x _i+1,j,c −x _i,j,c =0, and if j+1>N _h , x _i,j+1,c −x _{i,j, c} =0. Also, in equation (7), c is a value of 1, 2, or 3, and indicates red, green, or blue.

In equation (6), the function Ω _w (·) (where the subscript w is the weight vector w) is the ordered weighted L ₁ -norm (hereinafter “OWL” (Ordered Weighted L ₁ -norm)).

“Reference 1: X. Zeng and MA Figueiredo, “The ordered weighted L ₁ norm: Atomic formulation, projections, and algorithms”, arXiv:1409.4271, pp.1-13, 2014”

When a vector _{u =} ( ^u ₁ , _{u 2} , . expressed.

In equation (8), the operation |·| _[n] is an operation that selects the n-th largest value of |u ₁ |, |u ₂ |, ..., |u _N | when a vector u is given as an argument. is. In equation (8), the function |·| _↓ is given a vector u as an argument, and is a vector whose elements are the absolute values of the elements of the vector u, that is, (|u ₁ |, |u ₂ |, ..., | u _N |) is an operation that generates a vector in which the elements of ^T are sorted in descending order.

If the elements of the weight vector w of the function Ω _w (·) are non-negative and non _- decreasing, that is, (0≦w ₁ ≦w ₂ ≦ . become. In FIG. 1A, the contour line indicated by reference numeral 51 indicates the lowest position. The height of the contour line indicated by reference numeral 52 that is one outside the contour line indicated by reference numeral 51 is higher than the height of the contour line indicated by reference numeral 51. is higher than Hereinafter, an OWL with a non-negative, non-decreasing weight vector w will be referred to as a Pointed Star Norm (PSN) from the shape shown in FIG. 1(a).

The cusp norm is a non-convex function, as can be seen from the shape shown in FIG. 1(a). The cusp norm is a sparse regularization function with stronger sparsity regularizing properties than the L1 norm, which has the shape shown in Fig. ₁ (b), which is most commonly used in sparse regularization, by virtue of being a non-convex function. be. There is a norm called the _{Lp pseudo} -norm represented by the following equation (9) and having the shape of FIG. In FIGS. 1(b) and 1(c), as in FIG. 1(a), the contour lines on the outer side are higher than the contour lines on the inner side.

By using the cusp norm, the cusp gradient imposes a small weight regularization on pixels with large gradients. Therefore, in the cusp norm, a weight vector w represented by the following equation (10), which introduces two control parameters, a control parameter λ and a control parameter α, can be used.

In equation (10), vector 0 _N is a vector of length N with all elements "0" and vector 1 _N is a vector of length N with all elements "1". The first α number is 0, and the following N-α number is a vector of λ. A portion with a large difference among adjacent differences has a multiplier of 0 and is not included in the optimization target, and a portion with a small difference has a multiplier. becomes λ, which can be smoothed further. As a result, the smoothness of the gradation can be adjusted by the control parameter λ of the weight vector w, and the sharpness of the edge can be adjusted by the control parameter α. Note that the control parameter λ is a real number equal to or greater than 0, and the greater the value, the smoother the gradation can be adjusted. The control parameter α is an integer equal to or greater than 0, and can control the number of pixels with sharp edges.

In the weight vector of equation (10), when α=0, the weight vector w is obtained by multiplying the vector _1N by λ (w ₁ =w ₂ = . . . =w _N =λ). The optimization problem of equation ( ₃ ) with the cusp norm is equivalent to the L1 regularization problem formulated for TV.

On the other hand, in the case of (w ₁ =w ₂ =...=w _α =0, w _α+1 =w _α+2 =...=w _N =λ), if λ→∞, the equation (3 ) approaches the L ₀ pseudo-norm constraint problem formulated in the case of L ₀ gradient projection. Therefore, the cusp norm has the property of a sparse regularization function that is stronger than the L ₁ norm and milder than the L ₀ pseudo-norm.

For the control parameter λ that adjusts the smoothness of the gradation, it is necessary to adjust the value while viewing the output image in which the output image signal vector x is displayed on the screen. On the other hand, for the control parameter α for adjusting the sharpness of the edge, an integer of about several percent of the number N of pixels is determined.

As described above, since the cuspate norm is a non-convex function, it is known that finding the solution of the optimization problem of formula (3) generally requires a very large amount of computation time. In order to solve this computation time problem, in this embodiment, one of the algorithms included in the proximity separation method, the alternating direction multiplication method (hereinafter "ADMM" (Alternating Direction Method of Multipliers)).

"Reference 2: D. Gabay and B. Mercier, ``A dual algorithm for the solution of nonlinear variational problems via finite element approximation'', Computers and Mathematics with Applications, Vol.2 pp17-40, Pergamon Press, 1976"

In order to utilize ADMM, equation (3) needs to be transformed into a form in which ADMM can be applied. Therefore, a GPSL (Group Pointed Star L _1,2 -norm) function is defined as the following equation (11).

The vector y in equation (11) is a vector defined by the following equation (12).

In equation (11), the symbol represented by the following equation (13) is a symbol indicating a group for each pixel of n=1, . . . , N. Therefore, the argument of the function Ω _w (·) in equation (11) means that the vector y is divided pixel by pixel.

By using the GPSL function of equation (11), the function PSG _w (·) in the second term of equation (3) (where the subscript w is the weight vector w) is given by the following equation (14). can be transformed into

In Equation (14), matrix D is a matrix defined by Equation (15) below and is a discrete difference operator matrix with cyclic boundary conditions. In Equation (14), matrix B is a matrix defined by Equation (16) below, and is a 0-1 binary diagonal matrix with a cyclic boundary of 0 (see, for example, Non-Patent Document 3).

By the transformation shown in equation (14), the optimization problem represented by equation (3) can be transformed into the following equation (17) that can utilize ADMM.

By using the ADMM, it is possible to obtain the local optimum solution of the optimization problem of formula (17) based on the update formulas shown by the following formulas (18) to (20). In equations (18) and (19), γ is the step size and γ>0.

In the above equations (18) to (20), each of the equation (18), which is the update equation for the vector x ^(k+1) , and the equation (19), which is the update equation for the vector y ^(k+1) , is It's an optimization problem. Therefore, a solution method for obtaining each solution will be described below.

As for the vector x ^(k+1) , the following equation (21) is the solution as shown in Non-Patent Document 3.

In equation (21), the function F(·) is the two-dimensional discrete Fourier transform, and the function F ^* (·) is the two-dimensional inverse discrete Fourier transform. The matrix Λ is a diagonal matrix whose diagonal elements are the Fourier transform results of the Laplacian filter kernel. Therefore, the matrix of the following formula (22) included in the argument of F ^* (·) becomes a diagonal matrix, and the inverse matrix of formula (22) can be calculated by division for each element.

Next, the solution for vector y ^(k+1) will be described. First, the following formula (23) is defined.

Using the equation (23), the solution of the vector y ^(k+1) is expressed as the following equation (24).

The calculation of the following formula (25), which is the first term on the right side of the formula (24), can be expressed as the following formula (27) after replacing the following formula (26).

However, n=1, . . . , N in equation (27). In Equation (27), (p ₁ ^(k) , . . . , p _N ^(k) ) is expressed as Equation (28) below.

Let the argument of the function on the right side of the equation (28) be the vector v ^(k) as shown in the following equation (29).

The expression on the right side of expression (28) is defined as the following expression (30), and (p ₁ ^(k) , ..., p _N ^(k) ) is calculated by performing the operation on the right side of expression (30). can be done.

Here, since there is a relationship shown in the following equation (31), the calculation using the sign function sgn(·) in the right-hand side of equation (30) can be omitted.

The function represented by the following formula (32) included in the right side of the formula (30) indicates that the calculation of projection onto the non-increasing monotone cone represented by the following formula (33) is performed.

The operation of formula (32) is known as the Isotonic Regression problem. However, since the relationship shown in the following equation (34) exists, the projection by the equation (33) becomes an identity map and can be omitted.

As a result, formula (30) can be transformed into the following formula (35).

In equation (35), the function (·) ₊ is a ramp function that clips its argument to non-negative values. The matrix P ^(k) is a vector permutation matrix generated by sorting the elements of the vector v ^(k) in descending order, and has the relationship of the following equation (36).

As a result, by repeatedly executing the formula (21), the formula (24) applying the formulas (27) and (35), and the update formula of the formula (20), Can solve optimization problems. Note that the step size γ is updated as γ←ηγ (where 0<η<1) at each iteration. That is, updating is performed with the value obtained by multiplying the previous step size γ by the attenuation rate η as the new step size γ.

(Configuration of signal processing device)
FIG. 2 is a block diagram showing the internal configuration of the signal processing device 1 that calculates the solution of the optimization problem of equation (17). The signal processing device 1 includes a control unit 2. The control unit 2 includes a first variable update unit 10, a second variable update unit 20, a third variable update unit 30, an update determination unit 40, and a step size update unit. 50.

The update determination unit 40 stores the value of "k" (hereinafter referred to as "step k") indicating the number of times of update processing in an internal storage area. The update determination unit 40 outputs an instruction signal for updating the step size to the step size updating unit 50 at the timing of increasing the value of the step k. The update determination unit 40 determines whether or not to continue updating the vector x as the first variable, the vector y as the second variable, and the vector d as the third variable according to a predetermined end condition. do. Here, the termination condition is, for example, that the value of step k at that time coincides with a predetermined upper limit of the number of updates.

When the update determination unit 40 determines to continue the update process, it outputs the vector y ^(k) calculated one step before by the second variable update unit 20 to the first variable update unit 10, and performs the third variable update. The vector d ^(k) calculated by the unit 30 one step before is output to the first variable updating unit 10 , the second variable updating unit 20 and the third variable updating unit 30 . On the other hand, when the update determination unit 40 determines not to continue the process of updating, at that time, the vector x ^(k+1) calculated by the first variable update unit 10 is output to the outside as the output image signal vector x. output and terminate the process.

The step size updating unit 50 takes in the step size γ and the attenuation rate η given from the outside. The step size updating unit 50 writes and stores the acquired step size γ and attenuation factor η in its internal storage area. After writing the step size γ to the internal storage area, the step size updating unit 50 outputs the written step size γ to the first variable updating unit 10 and the second variable updating unit 20 . Upon receiving an instruction signal to update the step size from the update determination unit 40, the step size update unit 50 updates the step size γ stored in the internal storage area based on the attenuation rate η. The step size γ is output to the first variable updater 10 and the second variable updater 20 .

The first variable update unit 10 takes in an externally supplied input image signal vector f, a step size γ, a vector y ^(k) , and a vector d ^(k) . The first variable updating unit 10 calculates the vector x ^(k+1) , which is the first variable, by performing the calculation of Equation (21).

FIG. 3 is a block diagram showing the internal configuration of the first variable updating unit 10. As shown in FIG. The first variable update unit 10 includes a subtractor 101, a differential transpose calculator 102, an adder 103, a multiplier 104, an FFT (Fast Fourier Transform) unit 105, an inverse calculator 106, a coefficient generator 107, a multiplier 108 and an inverse An FFT unit 109 is provided.

Subtractor 101 subtracts vector d ^(k) calculated by third variable updating unit 30 from vector y ^(k) calculated by second variable updating unit 20 to calculate a subtraction value represented by the following equation (37). .

The difference transpose operation unit 102 stores the matrix D in advance in an internal storage area, and generates a matrix ^DT that is a transposed matrix of the matrix D stored in the internal storage area. The difference transpose calculation unit 102 multiplies the generated matrix ^DT and the subtraction value calculated by the subtractor 101 to calculate the multiplication value shown in the following equation (38).

The reciprocal calculator 106 takes in the step size γ and calculates the reciprocal of the step size γ, that is, γ ⁻¹ . Multiplier 104 multiplies the multiplied value calculated by differential transpose calculation section 102 and γ ⁻¹ calculated by reciprocal calculation section 106 to calculate the multiplied value. The adder 103 adds the input image signal vector f and the multiplied value calculated by the multiplier 104 to calculate the added value shown in Equation (39).

FFT section 105 performs two-dimensional fast Fourier transform on the added value calculated by adder 103 . Coefficient generation section 107 stores matrix I and matrix Λ in advance in an internal storage area, and based on matrix I and matrix Λ stored in the internal storage area and γ ⁻¹ calculated by reciprocal operation section 106 to generate the matrix of equation (22).

A multiplier 108 multiplies the output of the FFT section 105 and the matrix generated by the coefficient generating section 107 to calculate the multiplied value shown in the following equation (40).

The inverse FFT unit 109 performs two-dimensional inverse fast Fourier transform on the multiplied value calculated by the multiplier 108, that is, the calculation of equation (21) to calculate the vector x ^(k+1) as the first variable.

The second variable updating unit 20 updates the vector x (k+1), which is the first variable calculated by the first variable updating unit 10, and the vector x ^(k+1) , which is the third variable calculated one step before by the third variable updating unit 30. Take d ^(k) , the weight vector w, and the step size γ. The second variable update unit 20 calculates the vector y ^(k+1) , which is the second variable, by performing the calculation of the equation (24) applying the equations (27) and (35).

FIG. 4 is a block diagram showing the internal configuration of the second variable updating unit 20. As shown in FIG. The second variable update unit 20 includes a difference calculation unit 201, an adder 202, a binary mask unit 203, an L2 norm calculation unit 204, a descending order sort unit 205, a permutation matrix generation unit 206, a multiplier 207, a subtractor 208, and a ramp function unit. 209 , an inverse permutation unit 210 , a coefficient generation unit 211 , an element product calculation unit 212 , a binary mask unit 213 , a binary mask unit 214 and an adder 215 .

The difference calculation unit 201 stores the matrix D in advance in an internal storage area, and multiplies the matrix D stored in the internal storage area by the vector x ^(k+1) calculated by the first variable updating unit 10 to obtain Calculate the multiplication value. The adder 202 adds the multiplied value calculated by the difference calculator 201 and the vector d ^(k) calculated by the third variable updater 30 to calculate the added value shown in Equation (23). Binary mask section 203 stores matrix B in advance in an internal storage area, and multiplies matrix B stored in the internal storage area by the addition value calculated by adder 202 to obtain the result shown in equation (26). Calculate the multiplication value.

The L2 norm calculation unit 204 is a calculation unit that calculates the L2 norm for each pixel, and calculates the vector v ^(k) shown in Equation (29) from the multiplication value calculated by the binary mask unit 203. The descending order sorting unit 205 sorts the elements of the vector, whose elements are the absolute values of the elements included in the vector v ^(k) calculated by the L2 norm calculating unit 204, in descending order to obtain a vector a ^{( k)} .

Permutation matrix generation section 206 generates matrix P ^(k) which is a permutation matrix of vector a ^{( k) from vector v (k)} generated by L2 norm calculation section 204 . Multiplier 207 multiplies weight vector w and step size γ to calculate a multiplied value shown in the following equation (42).

Subtractor 208 calculates a subtraction value by subtracting the multiplication value calculated by multiplier 207 from vector a ^(k) generated by descending order sorting section 205 . The ramp function unit 209 performs an operation of clipping the subtraction value calculated by the subtractor 208 to a non-negative value, and calculates the operation result shown in the following equation (43).

Inverse permutation section 210 generates matrix P ^(k)T , which is a transposed matrix of matrix P ( ^k) generated by permutation matrix generation section 206 . The inverse permutation unit 210 multiplies the generated matrix P ^(k)T and the output of the ramp function unit 209 to calculate the vector p ^(k) given by the following equation (44).

^The coefficient generation unit 211 generates ^coefficients q _n ^{( k)} is calculated.

However, n=1, . . . , N in equation (45). The element product calculation unit 212 calculates a value represented by the following equation (46) based on the multiplied value represented by the equation (26) calculated by the binary mask unit 203 and the coefficient q _n ^(k) generated by the coefficient generation unit 211. Calculate the vector z ^(k)

However, in Equation (46), ⁿ =1, .

The binary mask unit 213 stores the matrix B in advance in the internal storage area, and multiplies the matrix B stored in the internal storage area by the vector z ^(k) calculated by the element product operation unit 212. Calculate the value. The binary mask unit 214 stores in advance the matrix (IB) obtained by subtracting the matrix B from the matrix I in the internal storage area. Calculates the multiplication value by multiplying the addition value calculated by . Adder 215 adds the output of binary mask section 213 and the multiplied value calculated by binary mask section 214 to calculate vector y ^(k+1) , which is the second variable shown in the following equation (48).

The third variable updating unit 30 updates the vector x ^(k+1) , which is the first variable calculated by the first variable updating unit 10, and the vector y ^(k+1) , which is the second variable calculated by the second variable updating unit 20. , and the vector d ^(k) calculated one step before by the third variable updating unit 30 . The third variable updating unit 30 calculates the vector d ^(k+1) , which is the third variable, by performing the calculation of Equation (20).

FIG. 5 is a block diagram showing the internal configuration of the third variable updating unit 30. As shown in FIG. The third variable updating unit 30 has a difference calculating unit 301 , an adder 302 and a subtractor 303 . The difference calculation unit 301 stores the matrix D in advance in an internal storage area, and multiplies the matrix D stored in the internal storage area by the vector x ^(k+1) calculated by the first variable updating unit 10 to obtain Calculate the multiplication value. The adder 302 adds the vector d ^(k) calculated one step before by the subtractor 303 and the multiplied value calculated by the difference calculation unit 301 to calculate an addition value. The subtractor 303 subtracts the vector y ^(k+1) calculated by the second variable updating unit 20 from the added value calculated by the adder 302 to calculate the vector d ^(k+1) which is the third variable.

(Processing by signal processing section)
FIG. 6 is a flow chart showing the flow of processing by the signal processing unit 1. As shown in FIG. The first variable update unit 10 and the update determination unit 40 of the signal processing unit 1 take in the input image signal vector f given from the outside. The second variable updating unit 20 takes in a weight vector w in which the externally provided control parameter λ and the control parameter α are predetermined. The step size updating unit 50 takes in the step size γ and the attenuation rate η given from the outside (step S1).

The update determination unit 40 multiplies the matrix D prestored in the internal storage area by the captured input image signal vector f to calculate the multiplied value shown in the following equation (49).

The update determination unit 40 sets the calculated multiplied values as the vector y ⁽⁰⁾ , which is the initial value of the second variable, and the vector d ⁽⁰⁾ , which is the initial value of the third variable. The update determination unit 40 initializes step k to "0" (step S2).

The update determination unit 40 determines whether or not the termination condition is satisfied (step S3). Whether or not the termination condition is satisfied is determined by whether or not the value of step k at that time coincides with a predetermined upper limit value for the number of updates. Here, it is assumed that the value of step k=0 does not match the predetermined upper limit of the number of updates. The update determination unit 40 determines that the end condition is not satisfied because the value of step k does not match the predetermined upper limit of the number of updates (step S3, No), and the initial value of the second variable A certain vector y ⁽⁰⁾ is output to the first variable updating unit 10, and the vector d ⁽⁰⁾ , which is the initial value of the third variable, is output to the first variable updating unit 10, the second variable updating unit 20, and the third variable It is output to the updating unit 30.

Hereinafter, for the sake of generalization, the value of the number of steps stored in the internal storage area by the update determination unit 40 is not limited to "k=0", but is assumed to be "k", and the subscript of each variable is (k ) or (k+1). In the first variable update unit 10, the subtractor 101 takes in the vector y ^(k) and the vector d ^(k) output by the update determination unit 40. FIG. Subtractor 101 subtracts vector d ^{( k) from vector y(k)} to calculate a subtraction value shown in equation (37). The subtractor 101 outputs the calculated subtraction value to the difference transposition calculation section 102 .

The differential transpose calculation unit 102 generates a matrix ^DT that is a transposed matrix of the matrix D stored in the internal storage area. The difference transpose calculation unit 102 multiplies the generated matrix ^DT and the subtraction value calculated by the subtractor 101 to calculate the multiplication value shown in Equation (38). The difference transposition calculation unit 102 outputs the calculated multiplication value to the multiplier 104 .

The reciprocal calculation unit 106, in parallel with the processing performed by the adder 101 and the differential transposition calculation unit 102, takes in the step size γ given from the outside. The reciprocal calculation unit 106 calculates γ ⁻¹ which is the reciprocal of the captured step size γ. Reciprocal calculation section 106 outputs the calculated γ ⁻¹ to multiplier 104 and coefficient generation section 107 . Multiplier 104 calculates the multiplied value by multiplying the multiplied value of equation (38) output from differential transpose calculation unit 102 and γ ⁻¹ output from reciprocal calculation unit 106, and the calculated multiplied value is added to the adder. 103. The adder 103 takes in an externally applied input image signal vector f. The adder 103 adds the fetched input image signal vector f and the multiplied value output from the multiplier 104 to calculate the added value shown in Equation (39). Adder 103 outputs the calculated addition value to FFT section 105 .

FFT section 105 performs a two-dimensional fast Fourier transform operation on the added value of equation (39) output from adder 103 and outputs the operation result to multiplier 108 . In parallel with the processing performed by the multiplier 104, the adder 103, and the FFT unit 105, the coefficient generation unit 107 generates γ ⁻¹ output from the reciprocal operation unit 106, the matrix I and the matrix Λ stored in the internal storage area. generates the matrix of equation (22). Coefficient generation section 107 outputs the generated matrix to multiplier 108 .

Multiplier 108 multiplies the calculation result output from FFT section 105 by the matrix of formula (22) output from coefficient generation section 107 to calculate the multiplied value shown in formula (40). Multiplier 108 outputs the calculated multiplied value to inverse FFT section 109 .

Inverse FFT section 109 performs two-dimensional inverse fast Fourier transform shown in equation (21) on the multiplied value of equation (40) output from multiplier 108 to calculate the first variable vector ^(k+1) . The inverse FFT unit 109 outputs the calculated vector ^(k+1) to the second variable update unit 20, the third variable update unit 30, and the update determination unit 40 (step S5).

In the second variable update unit 20, the difference calculation unit 201 multiplies the matrix D stored in the internal storage area by the vector x ^(k+1) output by the inverse FFT unit 109 of the first variable update unit 10 Calculate the value. Difference calculation section 201 outputs the calculated multiplication value to adder 202 . The adder 202 adds the vector d ^(k) output from the update determination unit 40 and the multiplied value output from the difference calculation unit 201 to calculate the added value shown in Equation (23). The adder 202 outputs the calculated added value to the binary mask section 203 and the binary mask section 214 (step S6).

Binary mask section 203 multiplies matrix B stored in an internal storage area by the added value output from adder 202 to calculate the multiplied value shown in equation (26). The binary mask unit 203 outputs the calculated multiplication value to the L2 norm calculation unit 204 and the element product calculation unit 212 (step S7). L2 norm calculation section 204 calculates vector v ^(k) shown in equation (29) from the multiplied value output from binary mask section 203 . L2 norm calculator 204 outputs the calculated vector v ^(k) to permutation matrix generator 206 , descending order sorter 205 , and coefficient generator 211 . The descending order sorting unit 205 sorts the elements of the vector, whose elements are the absolute values of the elements included in the vector v ^(k) output from the L2 norm calculating unit 204, in descending order to obtain the vector a ^(k) shown in Equation (41). to generate The descending order sorting unit 205 outputs the generated vector a ^(k) to the subtractor 208 (step S8).

The multiplier 207 takes in the weight vector w given from the outside and the step size γ in parallel with the processing of steps S6 to S8. Multiplier 207 multiplies weight vector w and step size γ to calculate a multiplied value shown in equation (42). The multiplier 207 outputs the calculated multiplied value to the subtractor 208 (step S9).

A subtractor 208 calculates a subtraction value by subtracting the multiplied value calculated by the multiplier 207 from the vector a ^(k) output by the descending order sorting unit 205 . Subtractor 208 outputs the calculated subtraction value to ramp function section 209 . The ramp function unit 209 performs the calculation shown in Equation (43) for clipping the subtracted value output from the subtractor 208 to a non-negative value, and outputs the calculation result to the inverse replacement unit 210 (step S10).

The permutation matrix generation unit 206 generates a matrix P ^(k) , which is a permutation matrix of the vector a ^{(k ) from the vector v(k)} output by the L2 norm calculation unit 204, in parallel with the process of step S10. The permutation matrix generation unit 206 outputs the generated matrix P ^(k) to the inverse permutation unit 210 (step S11).

Inverse permutation section 210 generates matrix P ^(k)T , which is a transposed matrix of matrix P ^(k) output by permutation matrix generation section 206 . The inverse permutation unit 210 multiplies the generated matrix P ^(k)T by the calculation result of the calculation of the formula (43) output by the ramp function unit 209, and obtains the multiplied value as the vector p ^{(k )} is calculated. The inverse permutation unit 210 outputs the calculated vector p ^(k) to the coefficient generation unit 211 (step S12).

Coefficient generating section 211 generates coefficient q _n ^(k) expressed by equation (45) based on vector v ^(k) output from L2 norm computing section 204 and vector p ^(k) output from inverse permuting section 210. Calculate The coefficient generator 211 outputs the calculated coefficient q _n ^(k) to the element product calculator 212 (step S13).

Element product calculation section 212 performs the calculation shown in formula (47) based on the multiplied value shown by formula (26) output from binary mask section 203 and the coefficient q _n ^(k) output from coefficient generation section 211. to calculate the vector z ^(k) shown in equation (46). The element product calculation unit 212 outputs the calculated vector z ^(k) to the binary mask unit 213 (step S14).

In parallel with the processing of steps S7 to S14, the binary mask unit 214 multiplies the matrix (IB) stored in the internal storage area by the addition value calculated by the adder 202 and represented by the equation (23). to calculate the multiplication value. Binary mask section 214 outputs the calculated multiplication value to adder 215 . In parallel with the processing performed by the binary mask unit 214, the binary mask unit 213 multiplies the matrix B stored in the internal storage area by the vector z ^(k) output by the element product operation unit 212, and obtains the multiplied value. calculate. Binary mask section 213 outputs the calculated multiplication value to adder 215 .

The adder 215 adds the multiplied value output from the binary masking unit 213 and the multiplied value output from the binary masking unit 214, and performs the operation shown in equation (48) to obtain the second variable vector y ^(k+1) . calculate. The adder 215 outputs the calculated vector y ^(k+1) to the third variable update unit 30 and the update determination unit 40 (step S15).

In the third variable update unit 30, the difference calculation unit 301 multiplies the matrix D stored in the internal storage area by the vector x ^(k+1) output by the inverse FFT unit 109 of the first variable update unit 10, and multiplies Calculate the value. Difference calculation section 301 outputs the calculated multiplication value to adder 302 . The adder 302 adds the vector d ^(k) output from the update determination unit 40 and the multiplication value output from the difference calculation unit 301 to calculate an addition value. Adder 302 outputs the calculated addition value to subtractor 303 . The subtractor 303 subtracts the vector y ^(k+1) output by the adder 215 of the second variable updating unit 20 from the added value output by the adder 302 to calculate the vector d ^(k+1) as the third variable. The subtractor 303 outputs the calculated vector d ^(k+1) to the update determination unit 40 (step S16).

The update determination unit 40 determines the vector x ^(k+1) output by the first variable update unit 10, the vector y ^(k+1) output by the second variable update unit 20, and the vector d ⁽ k+1) output by the third variable update unit 30. ^k+1) , the value of step k stored in the internal storage area is read out, "1" is added to the read value, and the value after the addition is set as the new value of step k, and stored in the internal storage area. Write down and memorize. The update determination section 40 outputs an instruction signal for updating the step size to the step size updating section 50 . Upon receiving an instruction signal to update the step size from the update determination unit 40, the step size update unit 50 multiplies the step size γ stored in the internal storage area by the attenuation factor η to calculate ηγ. The step size update unit 50 rewrites the step size γ stored in the internal storage area to ηγ in order to set the calculated ηγ as the new step size γ. After rewriting the step size γ in the internal storage area to ηγ, the step size updating unit 50 outputs the rewritten step size γ to the first variable updating unit 10 and the second variable updating unit 20 (step S17 ).

The update determination unit 40 performs the determination process of step S3 again. If the update determination unit 40 determines that the value of step k at that point in time matches the predetermined upper limit of the number of updates and that the end condition is satisfied (step S3, Yes), at that point , the vector x ^(k+1) taken in is output to the outside as the output image signal vector x (step S4), and the process is terminated.

On the other hand, if the update determination unit 40 determines that the value of step k at that time does not match the predetermined upper limit of the number of updates and does not satisfy the termination condition (step S3, No), At that point, the fetched vector y ^(k+1) is output to the first variable updater 10, and the vector d ^(k+1) is transferred to the first variable updater 10, the second variable updater 20, and the third variable updater. 30, and the processing after step S5 is performed again.

(Implementation results)
7 and 8 are diagrams showing the results of edge preserving smoothing performed using the signal processing apparatus 1 described above. FIG. 7 shows the result of implementation when the control parameter λ=50 for adjusting the smoothness of the gradation, and FIG. 8 shows the result of implementation when the control parameter λ=200. FIGS. 7(a) and 8(a) show results obtained when the control parameter α for adjusting edge steepness is set to approximately 8% of the number of pixels N, and FIGS. 7(b) and 8( b) is an implementation result when the control parameter α is set to about 4% of the number N of pixels.

In FIGS. 7 and 8, images indicated by reference numerals 60a, 60b, 70a, and 70b are images displayed on the screen of the output image signals obtained after performing the edge-preserving smoothing process under each condition. In FIGS. 7 and 8, the graphs indicated by reference numerals 62a, 62b, 72a, and 72b are respectively the line of reference numeral 61a in the image 60a, the line of reference numeral 61b in the image 60b, the line of reference numeral 71a in the image 70a, and the reference numeral 71b in the image 70b. is a graph obtained by plotting the values of the luminance signal on the line of . The solid lines in the graphs indicated by reference numerals 62a, 62b, 72a, and 72b indicate the luminance signal values of the input image signal, and the dotted lines indicate the luminance signal values of the output image signal. As can be seen from a comparison between FIGS. 7 and 8, smoothing is performed in which the smoothness of the gradation is lost as the control parameter λ increases. On the other hand, as can be seen by comparing FIGS. 7(a) and 7(b), or between FIGS. 8(a) and 8(b), as the control parameter α becomes smaller, the sharpness of the edge decreases. It can be seen that the smoothing performed by

The signal processing apparatus 1 according to the above embodiment adjusts the control parameter λ, which is a sparse regularization function of a non-convex function for the input image signal and is the first control parameter included in the weight of the non-convex function. is a sparse regularization function for a non-convex function that changes the smoothness of the gradation and changes the sharpness of the edge by adjusting the control parameter α, which is the second control parameter included in the weight of the non-convex function, by Apply to generate an edge-preserving smoothed output image signal from the input image signal by computing the solution of the optimization problem by global optimization. This allows both gradient smoothness and edge steepness to be adjusted independently in edge-preserving smoothing with global optimization.

The computational complexity of the flow chart shown in FIG. and the amount of calculation of O(N) required for the four arithmetic operations, and the amount of calculation for one iteration is O(NlogN). Therefore, the signal processing apparatus 1 of the present embodiment performs processing for solving an optimization problem including the _cuspate norm of a non-convex function. It is possible to do edge-preserving smoothing by global optimization with sparse regularization performance stronger than and milder than the L ₀ pseudo-norm, but as robust as the L ₁ norm . Here, having robustness equivalent to the L ₁ norm means that even if each element of the input takes a minute value compared to the L ₀ norm, it does not react sensitively and is appropriately sparse. It has the ability to regularize.

In the above-described embodiment, in the determination process of step S3 of FIG. 6, the update determination unit 40 sets the end condition that the value of step k at that time coincides with the predetermined upper limit of the number of updates. . On the other hand, the norm of the amount of change of the vector x ^(k) obtained by the following equation (50) may be equal to or less than a predetermined value as a termination condition for the determination processing in step S3 of FIG.

The following equation (51) is used as the objective function, and the change amount of the objective function obtained by the following equation (52) may be equal to or less than a predetermined value as the end condition of the determination process in step S3 of FIG.

In the above embodiment, the FFT unit 105 may perform a two-dimensional discrete Fourier transform instead of the two-dimensional fast Fourier transform, and the inverse FFT unit 109 may perform two-dimensional inverse fast Fourier transform, A two-dimensional inverse discrete Fourier transform may be performed.

In the above embodiment, the image signal to be processed is a color image signal having three channels of RGB. On the other hand, a color image signal having four channels of CYMK may be used as an image signal to be processed, or a monochrome image signal of one channel may be used as an image signal to be processed. In this case, in the case of RGB color image signals, N _c = ₃ ; in the case of CMYK color image signals, N _c =4; and in the case of monochrome image signals, N _c =1. Of the above equations, the following equations are modified as follows using the variable _Nc . The suffix “3N” in formulas (1) and (2) becomes “NN _c ”. The subscript “N _v ×N _h ×3” in Equation (5) becomes “N _v ×N _h ×N _c ”. Of the numerical values indicating the range of symbol Σ in equation (7), “3” above Σ becomes “N _c ”. The subscript "6N" in equation (12) becomes "2NN _c ". The subscript “6N×3N” in equation (15) becomes “2NN _c ×NN _c ”. The subscript “6N×6N” in equation (16) becomes “2NN _c ×2NN _c ”. In the definition formula that is the same as formula (2) under argmin of formula (18), the suffix "3N" becomes " _NNc " as in formula (2). In the definition formula that is the same as formula (12) under argmin of formula (19), the suffix "6N" becomes "2NN _c " as in formula (12).

In the above-described embodiment, the signal processing apparatus 1 performs edge-preserving smoothing on an image signal, that is, a two-dimensional signal, but a one-dimensional signal such as a communication signal may be processed. A similar effect can be obtained for

The signal processing device 1 in the above-described embodiment may be realized by a computer. In that case, a program for realizing this function may be recorded in a computer-readable recording medium, and the program recorded in this recording medium may be read into a computer system and executed. It should be noted that the "computer system" referred to here includes hardware such as an OS and peripheral devices. The term "computer-readable recording medium" refers to portable media such as flexible discs, magneto-optical discs, ROMs and CD-ROMs, and storage devices such as hard discs incorporated in computer systems. Furthermore, "computer-readable recording medium" means a medium that dynamically retains a program for a short period of time, like a communication line when transmitting a program via a network such as the Internet or a communication line such as a telephone line. It may also include something that holds the program for a certain period of time, such as a volatile memory inside a computer system that serves as a server or client in that case. Further, the program may be for realizing a part of the functions described above, or may be capable of realizing the functions described above in combination with a program already recorded in the computer system. It may be implemented using a programmable logic device such as an FPGA (Field Programmable Gate Array).

Although the embodiment of the present invention has been described in detail with reference to the drawings, the specific configuration is not limited to this embodiment, and includes design within the scope of the gist of the present invention.

It can be applied to techniques that perform edge-preserving smoothing through global optimization.

DESCRIPTION OF SYMBOLS 1... Signal processing apparatus 2... Control part 10... 1st variable update part 20... 2nd variable update part 30... 3rd variable update part 40... Update determination part 50... Step size update part

Claims (7)

1. By adjusting a first control parameter that is a sparse regularization function of a non-convex function for an input signal and is included in the weight of the non-convex function, the smoothness of the gradation is changed, and the weight of the non-convex function is adjusted. By adjusting a second control parameter contained in the input A signal processing apparatus comprising a controller for generating an edge-preserving smoothed output signal from a signal.
2. wherein the input signal and the output signal are image signals;
   The signal processing device according to claim 1.
3. The sparse regularization function of the non-convex function is the cusp norm defined by the following equation (1) and the elements of the weight vector w of the ordered weight L ₁ norm are non-negative non-decreasing (0 ≤ w1 ≤ w2 ... ≤ wN) is the cusp norm,
   

   

   3. The signal processing device according to claim 1 or 2.
4. The weight vector w is defined by the following formula (2),
   

   

   The first control parameter is λ in equation (2),
   The second control parameter is α in equation (2),
   The signal processing device according to any one of claims 1 to 3.
5. The optimization problem is
   Using the GPSL function defined by the following equation (3),
   

   

   defined as the following formula (4),
   

   

   The control unit calculates the solution of the optimization problem by an alternating direction multiplication method.
   5. The signal processing device according to any one of claims 1 to 4.
6. A signal processing method performed by a signal processing device,
   By adjusting a first control parameter that is a sparse regularization function of a non-convex function for an input signal and is included in the weight of the non-convex function, the smoothness of the gradation is changed and the weight of the non-convex function is adjusted. By adjusting a second control parameter contained in the input A signal processing method that produces an edge-preserving smoothed output signal from a signal.
7. A signal processing program for executing a computer as the signal processing device according to any one of claims 1 to 5.

Images