WO2001018743A1 - Fast and efficient computation of cubic-spline interpolation for data compression - Google Patents

Abstract

A fast and efficient method and system for computation of cubic-spline interpolation for data compression is described (II). In one aspect, the present invention is a method and system for obtaining a transform of the compressed signal and computing an inverse transform of the transform of the compressed signal to obtain the compressed signal (I).

Description

FAST AND EFFICIENT COMPUTATION OF CUBIC-SPLINE INTERPOLATION

FOR DATA COMPRESSION

FIELD OF THE INVENTION

This invention relates to data compression. More specifically, the invention relates to a new cubic-spline interpolation (CSI) for both 1-D and 2-D signals to sub-sample signal and image compression data.

BACKGROUND OF THE INVENTION

In most multimedia systems the amount of image data is so large that the use of image data compression is almost mandatory. Image data compression allows the image to be transmitted over the Internet in real time. Also it reduces the requirements for image storage. Presently both spatial and temporal data reduction techniques are available and continue to improve the performance of image data compression. The fundamental problem of image data compression is to increase the compression ratio and to reduce the computational complexity within an acceptable fidelity.

Interpolation is one of the more important functions that can be used in the process of estimating the intermediate values of a set of discrete sampling points. Interpolation is used extensively in image data compression to magnify or reduce images and to correct spatial distortions. For example, see R. G. Keys, "Cubic Convolution Interpolation for Digital Image Processing," IEEE Trans, on Acoustics, Speech, and Signal Processing, vol. ASSP-29, no.6, pp.1 153-1 160, Dec. 1981 , [1], the contents of which are hereby expressly incorporated by reference. In general, the process of decreasing the data rate is called decimation and the process of increasing data samples is called interpolation as described in H. S. Hou, and H. C. Andrews, "Cubic Splines for Image Interpolation and Digital Filtering," IEEE Trans, on Acoustics, Speech, and Signal Processing, vol. ASSP-26, no.6, pp.508-517, Dec. 1978, [2], the contents of which are hereby expressly incorporated by reference. It is well known that several interpolation functions such as linear interpolation (see,

W. K. Pratt, Digital Image Processing, second edition, John Wiley & Sons, Inc., New York, 1991 , [3], the contents of which are hereby expressly incorporated by reference.) cubic-convolution interpolation (see [1], and [3]), cubic B-spline interpolation (described in C. de Boor, A Practical Guide to Splines. New York: Springer-Verlag, 1978, [4]; M. Unser, A. Aldroubi, and M. Eden, "B-Spline Signal Processing: Part II-Efficient Design and Applications," IEEE Trans, on Signal Processing, vol.41, pp.834-848, Feb. 1993, [5]; M. Unser, A. Aldroubi, and M. Eden, "Enlargement or Reduction of Digital Images with Minimum Loss of Information," IEEE Trans, on Image Processing, vol.4, pp.247-258, March 1995, [6]; and [2]) can be used in the image data compression process.

The disadvantage of these interpolation schemes is that in general they are not designed to minimize the error between the original image and its reconstructed image. In 1981 Reed (I. S. Reed, Notes on Image Data Compression Using Linear Spline

Interpolation, Department of Electrical Engineering, University of Southern California, Los Angeles, California, 90089-2565, U.S.A., Nov. 1981 [7], the contents of which are hereby incoφorated by reference) and in 1998 Reed and Yu (I. S. Reed and A. Yu, Optimal Spline Interpolation for Image Compression, United States Patent, No. 5822456, Oct. 13, 1998 [8], the contents of which are hereby incoφorated by reference) developed a linear spline inteφolation scheme for re-sampling the image data. This linear spline inteφolation is based on the least-squares method with the linear inteφolation function.

Using an extension of the ideas of Reed in [7,8], a modified linear spline inteφolation algorithm, called the cubic-spline inteφolation (CSI) algorithm, is developed in this invention for the sub-sampling of image data. (The linear spline inteφolation explained in [8] and used by America On Line™ (AOL) will be called the "AOL algorithm" in this document from hereon.)

It follows from [1] that the cubic-convolution inteφolation, which is different from the B-spline inteφolation, can be performed much more efficiently than that of the cubic B-spline inteφolation method. In this invention, the new CSI scheme combines the least- squares method with a cubic-spline function developed by Keys [1] for the decimation process. Also the cubic-spline reconstruction is used in the inteφolation process. Therefore, the CSI constitutes a new scheme that is quite different from both cubic B- spline inteφolation [2,4-6] and cubic-convolution inteφolation [1,3]. The concept of the CSI for both 1-D and 2-D signals is describes and demonstrated in the following sections. In addition, it is shown by computer simulation that the CSI scheme obtains a better subjective quality for the reconstructed image than linear inteφolation, cubic-convolution inteφolation, cubic B-spline inteφolation and linear spline inteφolation. An important advantage of this new CSI scheme is that it can be computed by a use of the FFT technique. The complexity of the calculation of the CSI scheme is substantially less than other conventional means.

W. B. Pennebaker and J. L. Mitchell, JPEG Still Image Data Compression Standard, Van Nostrand Reinhold, New York, 1993, [9], the contents of which are hereby incoφorated by reference, describes the JPEG still image data compression standard. It is well known that the JPEG (see [9]) algorithm is the international compression standard for still-images. The disadvantage of the conventional JPEG algorithm is that it causes visually disturbing blocking effects when high quantization parameter is used to obtain a high compression ratio. One embodiment of this invention includes a simpler and modified JPEG encoder-decoder to improve the JPEG standard with a high compression ratio and still maintain a good quality reconstructed image.

Recently, the authors in T. K. Truong, L. J. Wang, I. S. Reed, W. S. Hsieh, and T.C. Cheng "Image data compression using cubic convolution spline inteφolation," accepted for publication in IEEE Transactions on Image Processing [10], the contents of which are hereby incoφorated by reference, proposed the modified JPEG encoder-decoder for τ= 2 that utilizes the CSI scheme with a compression ratio of 4 to 1 as a pre-processing stage of the JPEG encoder and the cubic-spline reconstruction with a ratio of 1 to 4 as a post- processing stage of the inverse JPEG decoder to achieve a high compression ratio.

In such a modified JPEG encoder the CSI scheme is the pre-processing stage of the JPEG encoder. It can be implemented by the use of the FFT algorithm. In addition, the output of the modified JPEG encoder represents the compressed data to be transmitted. It can be pre-computed and stored. In such a modified JPEG decoder, the cubic-spline reconstruction constitutes the post-processing stage of the JPEG decoder. This postprocessing stage is different from the conventional post-processing algorithms that were proposed to reduce the blocking effects of block-based coding in B. Ramamurthi and A. Gersho, "Nonlinear space variant post-processing of block coded images," IEEE Trans, on Acoustics, Speech, Signal Processing, vol. ASSP-34, pp.1258-1267, 1986, [11], Y. Yang, N. Galatsanos, and A. Katsaggelos, "Projection-based spatially adaptive reconstruction of block-transform compressed images," IEEE Trans, on Image Processing, vol.4, pp.896-908, July 1995 [12], the contents of which are hereby incoφorated by reference.

The proposed post-processing stage is an inteφolation process that uses the cubic- convolution inteφolation. In [10], a computer simulation shows that the modified JPEG encoder-decoder for τ= 2 obtains a better subjective quality and an objective PSNR of the reconstructed image than the JPEG algorithm described in T. Lane, Independent JPEG Group's free JPEG software, 1998, [13], the contents of which are hereby incoφorated by reference; and [9]. Furthermore, the modified inverse JPEG decoder requires less computational time than the conventional JPEG decoder. But, the disadvantage of the modified JPEG encoder-decoder for τ= 2 is that the computational time required for the modified JPEG encoder is greater than the conventional JPEG encoder.

Thus, in one aspect, the present invention describes a fast method to compute the modified JPEG encoder. It is shown in this aspect of the invention that the speed of the new method for computing the modified JPEG encoder is approximately two times faster than that of the conventional JPEG encoder with still a good quality of reconstructed image. SUMMARY OF THE INVENTION

The present invention describes a new cubic-spline inteφolation (CSI) for both 1-D and 2-D signals to sub-sample signal and image compression data. This new inteφolation scheme which is based on the least-squares method with a cubic-spline function can be implemented by the fast Fourier transform (FFT), and/or by a Winograd discrete Fourier transform (WDFT). The result is a simpler and faster inteφolation design than can be obtained by conventional means. It is shown by computer simulation that such a new CSI yields the most accurate algorithms for smoothing. Linear inteφolation, linear spline inteφolation, cubic-convolution inteφolation and cubic B-spline inteφolation tend to be inferior in performance. In addition it is shown in this invention that the CSI scheme can be performed by a fast and efficient computation. The proposed method uses a simpler technique in the decimation process. It requires substantially fewer additions and multiplications than the original CSI algorithm. In one aspect, the present invention is a method and system for defining a cubic- spline filter; correlating the filter with the signal to obtain a correlated signal; autocorrelating the filter to obtain autocorrelated filter coefficients; computing a transform of the correlated signal and the autocorrelated filter coefficients; dividing the transform of the correlated signal by the transform of the autocorrelated filter coefficients to obtain a transform of a compressed signal; and computing an inverse transform of the transform of the compressed signal to obtain the compressed signal. The signal, the filter, and the transforms may be one dimensional or two dimensional. Further, the transforms may be a fast Fourier transform (FFT) or a Winograd discrete Fourier transform (WDFT) with an overlap-save scheme. Also, a zonal filter may be defined to simplify the steps of correlating and autocorrelating.

Furthermore, a new type of overlap-save scheme can be utilized to solve the boundary-condition problems that occur between two neighboring sub-images in the actual image for higher compression ratios. It is also shown in this invention that a very efficient 9-point Winograd discrete Fourier transform (WDFT) can be used to replace the FFT needed to implement the CSI scheme image for higher compression ratio of 9 to 1. Finally, a fast new CSI algorithm is used along with the Joint Photographic Experts Group (JPEG) standard to design a modified JPEG encoder-decoder for image data compression. As a consequence, for the higher compression ratios the proposed modified JPEG encoder-decoder obtains a better quality of reconstructed image and also requires less computational time than both the conventional JPEG method and the America on Line (AOL) algorithm. BRIEF DESCRIPTION OF THE DRAWINGS

The objects, advantages and features of this invention will become more apparent from a consideration of the following detailed description and the drawings, in which:

FIG. 1 is an exemplary periodic function with a period of τ= 5;

FIG. 2 is an exemplary 1-D cubic-spline function;

FIG. 3 is an exemplary shifted cubic-spline function for n=6;

FIG. 4 is an exemplary side view of a 2-D cubic-spline function; FIG. 5 is an exemplary reconstructed function between sampling periods;

FIGs. 6 (a) - (d) are exemplary zonal masks of the 2-D cubic-spline function used in the computation of Y _Mϊ in (25) and A_jχh^_kι in (26) for the CSI scheme;

FIG. 7 is an exemplary overlap map grid points for calculating c. for τ=3; FIG. 8 is an exemplary overlap map grid points for calculating β for τ=3; FIG. 9 is an exemplary overlap map grid points for calculating γ for τ=3;

FIG. 10 is an exemplary overlap map grid points for calculating δ for τ=3; FIG. 1 1 is a 19^χ 19 sub-images of size 9><9 in an exemplary image of size 171 ^χ171; FIG. 12 is a reconstructed image with serious artifacts that is generated by using FCSI method implemented by the direct use of the 9 * 9 Winograd DFT for compression; FIG. 13 is an illustrative example of the FCSI algorithm implemented by a 5^χ5

Winograd DFT and the overlap-save sub-image method;

FIG. 14 is a reconstructed image with no apparent artifacts that is generated by using a FCSI implemented by the 9 <9 Winograd DFT and the overlap-save method for compression; FIG. 15 is an exemplary modified JPEG encoder for τ= 2, 3

FIG. 16 is an exemplary modified JPEG decoder for τ= 2, 3;

FIGs. 17 (a) - (d) illustrate some reconstructed images with a compression ratio of 100: 1 ; and

FIGs. 18 (a) - (d) illustrate some reconstructed images with a compression ratio of 200:1.

DETAILED DESCRIPTION

The present invention describes a new cubic-spline inteφolation (CSI) for both 1-D and 2-D signals to sub-sample signal and image compression data. This new inteφolation scheme which is based on the least-squares method with a cubic-spline function can be implemented by the fast Fourier transform (FFT). Furthermore, in one embodiment of the present invention, a new type of overlap-save scheme is utilized to solve the boundary- condition problems that occur between two neighboring sub-images in the actual image for higher compression ratios. In one embodiment of this invention, an efficient 9-point

Winograd discrete Fourier transform (WDFT) is used to replace the FFT needed to implement the CSI scheme image for higher compression ratio of 9 to 1. Finally, a fast new CSI algorithm is used along with the Joint Photographic Experts Group (JPEG) standard to design a modified JPEG encoder-decoder for image data compression.

In order to speed up the modified JPEG encoder for τ- 2 with a good image quality, the compression ratio of the CSI scheme is extended to 9 to 1 (τ= 3). However, if the compression ratio for the CSI scheme is increased, there is need for additional computations that involve considerably more additions and multiplications. Therefore, in this invention a new faster and efficient algorithm for CSI is developed. The basic idea of this new algorithm, called the fast cubic-spline inteφolation (FCSI) scheme, is based on the CSI scheme, but it has a simpler form than that used for the original CSI scheme. The FCSI scheme substantially reduces the complexity of the additional computations that are required for the increased compression ratio. Moreover, the constants a,β,γ and δ , used to compute the terms A _k^k for the new FCSI scheme, are accurately calculated in detail in this invention.

For τ= 3, the actual image to be compressed, a novel overlap-save sub-image method is used to solve for the needed boundary conditions. Also a simplified and efficient algorithm that uses a 9-point Winograd DFT (WDFT) instead of the FFT is implemented to compress an actual image. Computer runs show that for some gray images of size 512 by 512, the computational time of the FCSI encoder that is implemented by a new type of overlap-save sub-image method and the 9-point Winograd DFT on a 400-MHz Intel Pentium II personal computer using C-code can be dramatically reduced. The FCSI encoder requires only around 0.15 sec when compared with around 0.57 sec for the original CSI encoder at a compression ratio of 9 to 1 (τ= 3). Also, the new FCSI scheme obtains a PSNR that is similar to all of the other more complicated zonal filters considered in this invention. Finally, this FCSI scheme with τ= 3 is combined with the JPEG standard to speed up the modified JPEG encoder for color image coding and still obtain a better quality of reconstructed image than the JPEG algorithm for higher compression ratios. In other words, the modified JPEG encoder for τ= 3 requires 0.71 sec, 0.38 sec and 0.67 sec less time than the modified JPEG encoder for τ= 2 [10], the conventional JPEG encoder [13] and the AOL algorithm [8], respectively.

This document is organized as follows: In Section II the encoding algorithm for this new inteφolation method is derived in detail. In addition it is shown now that the properties of the FFT and the convolution theorem can be used to compute the CSI scheme. The decoding algorithm is explained in Section III. In Section IV the FCSI computation is described. The constants needed for the FCSI are calculated in Section V. A novel FCSI algorithm that uses a new type of overlap-save sub-image technique and the

Winograd DFT algorithm is developed in Section VI. In Section VII the modified JPEG encoder-decoder is presented. Finally, experimental results are presented in Section VIII. II. ENCODING ALGORITHM FOR THE NEW INTERPOLATION METHOD

Encoding with the CSI scheme, utilizes the decimation process needed to perform image data compression. The philosophy of the CSI scheme is to recalculate the sampled values of the signal or image data by means of the least-squares method using the cubic- spline function. It is shown in this section that this new proposed method applies to both 1-D and 2-D signals as follows: A. CSI for the 1-D Signal

Let τ be a fixed, positive integer. Also let the data function X(t) be periodic with period nτ , where n is an integer. An example of X(t) for n = 5 is illustrated in FIG.1 From [1], the 1-D cubic-spline function R(t), shown in FIG.2, is defined by l

R(t) 2 (1)

Next one needs to define the shift function of the cubic-spline function R(/) as follows: Ψ_k(t) = R(t-kτ) for 0 ≤k≤n -1. (2)

An example of the shift functions Ψ_k (t) for n = 6 is shown in FIG.3. The goal is to approximate X(t) by a .-point sum, given by

S(t) = ∑X_kV_k(t) = ∑X_kR(t-kτ) (3) k=0 k=0 in a least-squares fashion, where X₀,---,X_n_^ are the reconstructed values at the sampling points which represent the compressed data to be transmitted or stored. The function S(t) in (3) is cubic-spline reconstruction of the function X(t) using the weights X₀ , ^{■ ■ ■} , X_n__x . It follows from [7] that the least-squares approximation of S(/) to X(t) is defined by

L(X₀ ,X_lt- X,^ ) = ∑ (X(t) - 5(t))² = ∑ {X( - ∑ X_k Ψ_k (Of , (4) r=-2r l=-2τ *=0 where the sum is taken over one period nτ plus 2r of the data. FIG. 1 shows an exemplary periodic function with a period 5r, FIG.2 is an exemplary 1-D cubic-spline function, and FIG.3 is an exemplary shifted cubic-spline function for n = 6.

Using the same procedure, described in [7,8], one can find the weights X₀,X_ly---,X_n„_s in (4) that minimize the function L(X₀,X_l,---X_n__l). To minimize (4) partial differentiation of L(X₀,X_l,---X„__l) with respect to Jζ for 0<J≤n-\ yields the following set of equations: c%(X₀,X -X^)

^■= ∑ 2(X(t)-∑X_kV_k(t))-Ψ_j(t) = 0 dX. k=0

-1- or

∑A _kX_k =Y, for O≤J≤n-l , (5) where

'/.A ∑Ψ_λ(t)Ψ,( f r 0≤j,k≤n-l, (6) and

Y_J = ∑X(t)Ψ_t(t) for O≤ ≤n-1. (7)

.=-2r

The term Y in (7) can be reduced as follows: Y,= ∑ X(t)R(t - jτ) = ^U∑X(t)R(t - jτ) . (8) i=-2r '=0-2)r

Let t - jτ = m, then

Y,= ∑X(m ₊ jτ)R(m). (9) m=-2r+l

Note that computing Y ^' in (9) involves n correlation coefficients of only 3r-l points.

Now let the periodicity of the cubic-spline function beR(t) = R(t + nτ) , i.e., R(t) has a period of nτ. Then the matrix form of A _k in (6) can be simplified to: J.k ∑' R(m + (j- k)r)R(m)

_∑2r- !_2τ+i(R(m))- =« if Omodn

∑^2r- , R(m + τ)R(m) = β if j - k ≡ ±1 mod n (10) lr-_2τ+lR(m + 2τ)R(m) = y if j - k ≡ ±2 mod n ∑J^~l_2τ+, R(m + 3τ)R(m) = δ if j - k ≡ ±3 mod n 0 otherwise

It follows from [7] that A _k in (10) can be expressed in circulant form by letting

^AJ,k ^{~ B}(k-_J)„ i (11) where (k - j)_n denotes the residue (k - j) modulo. , and

5₀=α,5,= ?,5₂ =_r,B₃=δ,B₄=0, -, B_n_₄=0,B_n_₃ =δ,B„_₂=γ, (12)

B.. β-

Therefore, A _k in (11) and (12) has the following symmetric, circulant representation:

A substitution of (13) into (5) yields the matrix equation,

A-X = Y, (14) where the matrix A is given in (13), X = (X_Q , X_χ , ^■ • • , X_n_ f and Y = (Y₀,Y, — ,Y,,__l)^τ .

Since the /. χ« matrix on the left side of (14) is a circulant matrix, (14) reduces immediately to

* =∑ ,>. for0<y<«-l. (15)

But B_t =5„_, for = 1,2, - • • , j_« / 2 J , where \_χ\ denotes the greatest integer less than or equal to x. Thus (15) becomes r,=∑ „*.,-*>, for O≤j≤n-l. (16)

In (16) the FFT can be used to solve for the X_k . To see this, let the FFT of Y_} , X_k and B_j for 0<j,k,m≤n-\ be defined by Y_m , X_m and B_m , respectively. By using the convolution theorem describes in E. O. Brigham, The Fast Fourier Transform and its Application, Prentice-Hall International, Inc., Englewood Cliffs, New Jersey, 1988 [14]; and A. V. Oppenheim, and R. W. Schafer, Digital Signal Processing, Prentice-Hall, Inc., Englewood Cliffs, New Jersey, 1975 [15], the contents of which are hereby expressly incoφorated be reference, it is easy to see that the solution of (16) in the frequency domain can be expressed as Y_m = X_m ^■ B_m or X_m = Y B_m , where B_m ≠ 0. Thus, using the inverse FFT of X_m , the X_k can be obtained for 0 < k < n - 1.

The encoding method for the 1-D signal is summarized as follows:

Choose a suitable value of integer r . The compression ratio is roughly r . Apply (9) to find Y_} . Also apply (10), (11) and (12) to find B_} .

Find the FFT of Y_t and _3_y to obtain Y_m and B_m, respectively. Also compute

^Xm ~^Ym^IBm-

Take the inverse FFT of X_m to get X_k which is the compressed data to be transmitted or stored. 1 B. CSI for the 2-D Signal

Let X(t_l,t₂) be a doubly periodic signal (e.g., image) of periods n_xτ and n₂τ with respect to integer variables t, and t₂ , where /., and n₂ are also integers. A 2-D cubic- 5 spline function, R(t, , t₂ ) , is defined by

R(t_l,t₂) = R(t_l)-R(t₂), (17) where R(.,) and R(t₂) are 1-D cubic-spline functions, respectively. A 3-D plot of this cubic-spline function is shown in FIG.4. It is well known fact that 2-D inteφolation can be accomplished by the use of 1-D inteφolations with respect to each coordinate [1,3]. 10 By analogy with the 1-D case, let

= R(t_x-k_xτ)-R(t₂-k₂τ) for 0 ≤ k, ≤ n, -1 and / = 1,2. (18)

By a procedure similar to that used in the 1-D case in (3), one defines the 2-D CSI as follows:

15

S(tf₂)=∑ 5' Ψ , (',,',) > (19)

A,=0 *₂=0 where X_{k Aι} are the reconstructed values at sampling points which represent the compressed image to be transmitted or stored. Again one wants to find the best weights X_{k k} such that

20 W ,)= Σ Σ (^ '2)- ∑ ∑ Ψ ,(',,'₂))² (20)

,,=-2r .₂=-2r i_|=0 *,=0 is a minimum. It follows from [7] that minimizing (20) yields

~ ' ² - ∑ Σ 2(X(t,t₂)-∑ ∑^ Ψ ,(/„/₂))Ψ_ΛJl(/₁,/₂) = 0,

/,=-2r .,=-2r i,=0 *₂=0

25 or

A- Σ,=0 * Σ,=0

^{> 0≤}J. ^≤n,^{~l and = 1}'²- (²¹) where ,_., _ = ∑ Σ ^Ψ _ΛΛ('..'2)Ψ _.(',.'2). 0≤ _ι,Λ,<n_I-l and = l,2, (22)

.,=-2r r,=-2r

30 and

> _.2 = Σ Σ ('. ^2 )Ψ,.„ Ci ,'2 ). 0 < < «, - 1 and i = 1, 2. (23)

.,=-2. ,=-2r

The term y in (23) reduces as follows:

Let / _A - j_kτ = m_k for / = 1,2 , then 2r-l ^γ _{h j}, = Σ Σ X^'(m_l+j,τ,m₂+j₂τ)R(m_l, ₂), (25) m,=~2τ+\ m₂=-2r+l where R(/.z,,/..₂) is the 2-D cubic-spline function shown in FIG.4. The terms A _k^ in (22) are treated in a similar fashion as follows:

= Σ Σ R(m_i+(j_]-k_])τ,m₂+(j₂-k₂)τ)R(m_l,m₂) (26)

if , -A, ≡0mod«_| andy₂ -k₂ ≡0mod«₂

:τ,m₂)R(m ,m₂) = β ≡±lmodn, and_₂ -A, ≡0 odn ≡0mod«, andy₂ -k₂ ≡±lmodn,

2-,/,,. »,,=-2.+l "''' ± 2r, JH . )R(m_l , m₂ ) = / ≡±2mod«_! andy₂ -k₂ ≡0mod .₂ ≡Omodn, and_/₂ -k₂ ≡±2modn₂

/.(/, ± 3r, «₂ )R(m_l , .«₂ ) = δ ≡ ±3 mod/?, and j₂ -k₂ ≡0moάn₂ ≡ Omodπ, andy₂ ~k₂ ≡ ±3mod/?₂

Σ^ , (m ± τ,n ±τ)R(m_x , //?, ) = η ≡±lmodϊ, and^, -k₂ ≡±lmod//, 2r-I

∑² -_r ₊ l ^R('ⁿ ± 2r, m₂ ± r)R(m_t , m₂ ) = λ ≡±2mod«, andy, -k₂ ≡±lmod /?₂ ≡ ±lmod /_| andy₂ -k₂ ≡ ±2 mod ?₂

±3r,m₂ ± τ)R{m_x , w?₂ ) = ≡ ±3 mod/?, andy₂ -k₂ ≡±lmod.?, ≡±lmod/z, andy₂ ~k₂ ≡±3mod/?₂

^■-2r-l

, Λ(« ±2τ,m₂±2τ)R(m_l,m₂) = p ≡ ±2 mod «, and j₂-k₂ ≡ ±2 mod /?₂ 2r-l

-2r+l /.(//?, ±±33τr,,mw₂₂ ±±22ττ))RR((mm_i,mm₂₂)) == σσ

≡ ±3 mod ι, and j₂ - A-, ≡±2 mod n₂ oτj_i -A, ≡±2mod/?₁ and j₂ -A, ≡±3mod/?₂ 0 otherwise

From Eq. (26) given above, the array A _it can be expressed in 2-D circulant form as follows:

^AJU2 ^2 ~^B ^-ιύ,_h, (*2-/2.„₂ ' (27) where (A, -_/,)„ denotes the residue (A, -7,) modulo /., for . = 1,2 and

I l

O 01/18743

where 0 <_.,</.,- 1 for / = 1,2. Note that if the array [_5_ ] is represented in matrix form, it is a block-circulant matrix. FIG.4 is the side view of the 2-D cubic-spline function. Since the matrix [B_s ] in (28) is a block-circulant matrix, (21) can be expressed by

y_h._h = ∑'∑'^Λ,-_.,),,, a_.-_..)._,. 0≤Λ ≤«, -1 and / = 1,2, (29) where (A, -7,),, denotes (A, -y.)mod/ι, for . = 1,2. Using a procedure similar to that used in (16), then (29) becomes

I ,. _-A.>„_.' O≤ ≤".-¹ and ¹ = 1,2. (30)

In (30) the 2-D FFT [14,15] can be used to solve for the X_k Let the 2-D FFT of

Λ-/_ ' _Y _ιand B_Sι for

be defined, respectively, by Ϋ_m „ , _?,„ „ and _5_m „ for 0 </«<«, -1 , 0 < n < n₂ -1. Then Eq. (30) in the frequency domain can be expressed as Y_m,„ =X,„_,„ -B_mn or X_m,n=Ϋ_m, B_mn, where _?„,,, ≠O. Finally, using the 2-D function of the inverse FFT of X_m<„, the X_{k k} are obtained for 0 < A, < /., - 1 , = 1,2.

The encoding method for the 2-D signal is summary in the following steps:

• Choose a suitable value of integer τ . The compression ratio is roughly r².

• Apply (25) to find Y_ΛJι . Next, apply (26), (27) and (28) to find B,_tJι .

• Take the 2-D FFT of Y and _5 to obtain Y_mn and B_mn, respectively. Also compute X_mn =Y_MIB_mn.

• Take the inverse FFT of X _mn to obtain X _{k k} which is the compressed image to be transmitted or stored.

III. DECODING ALGORITHM

In the decoding process, using the reconstructed values at the sampling points (e.g., X _k and X_{k Λ} ) obtained in Section II, the reconstructed points between the sampling points are obtained by means of the cubic-spline function. This decoding algorithm is called the cubic-spline reconstruction.

A. Decoding of the compressed 1-D Signal Since the n reconstruction values X₀ , ^{■ ■ ■} , X_n_ are known, the reconstructed signal

S(.) can be obtained by the use of (3). In other words, the retrieved signal is the convolution of the cubic-spline function R(/), defined in (1), and the sequence of n reconstructed values with sampling interval τ . The reconstructed function S(t_α ) between the two adjacent reconstructed values X_k and X_k+1 is illustrated in FIG.5 and given by the sum,

S(t_a) = X_k__xV_k__x(t_a) ₊ X_kV_k(t_a) ₊ X_k+V_k+x(t_a) ₊ X_k+2V_k+2(t_a) , (31) where kτ <.. <(A + l)r and Ψ_k(t) is defined in (2) and the boundary conditions as given in [1] are X_ =3(X₀ -X,) + X₂ and X„ =3(X_n_ - X_n_₂) +X„_₃. FIG. 5 is a reconstructed function between the sampling points.

B. Decoding of the compressed 2-D Signal

Since the reconstruction values X_{k k} for 0</, ≤n, -1,.=1,2 are known, the 2-D reconstructed image S(/,,t₂) can be obtained by the use of (19). In other words, the retrieved image is the 2-D convolution of the 2-D cubic-spline function R(/,,.,), given in

(17), and the 2-D sampled waveform X_{k k} . It is well known that the computationally simpler method described in [3], called bilinear inteφolation, can be used also to perform the 2-D inteφolation. The ideas of the bilinear inteφolation can be used to accomplish the 2-D cubic-spline reconstruction. In other words, the discrete data of each row can be inteφolated from the reconstructed values X_{k Λ} with similar inteφolation for the given discrete data of each column. It is easy to show that the reconstructed image S(t, ,/ ) between the four adjacent reconstructed values is given by

S(/,,/,) = ∑ ix_kl+mk2+nV_kl+mk2+n(t,dj) > (32) m=-\ n=-\ where A, r < /, < (A, + l)r , A, r < 1₁ <(k₂ + l)r and the boundary conditions as given in [1] are

X n-2.n) ⁺ X n-l,n ^•

IV. FAST COMPUTATION OF CSI

In fact, the CSI scheme needs a large number of pixels of the 2-D cubic-spline function R(t,,t₂) in order to compute the Y _hh in (25) and the A_{hJ kk} in (26). For example, computing each Y_Mι in (25) and A_hJι<kkι in (26) for τ= 2 (compression ratio 4: 1) and τ= 3 (compression ratio 9: 1) involve 81 and 169 pixels of the domain of function Λ(.,,.₂), respectively. As a consequence, the complexity of the calculation of (25) and (26) needed in the CSI scheme is substantially increased when the compression ratio is extended from τ= 2 to τ= 3. To illustrate this, for τ = 3 some abbreviated solutions are used to overcome the above problem of extra computational complexity. First, one defines a zonal mask of 169 pixels of ?(., , ., ) shown in FIG. 6(a), called zonal filter 1. This means that to compute each Y in (25) one needs to use a zonal mask of 169 pixels of R(t, ,t₂) to correlate with a periodic image, X(t_] ,t₂ ) . Also, to calculate each A _kιk in (26) one needs to use the 169 pixels for the function R(t, ,t₂) to autocorrelate with the other 169 pixels of R{t ,t₂) .

FIGs. 6(a) -(d) are zonal masks of the 2-D cubic-spline function used in the computation of Y in (25) and A_{JJι <k k^} in (26) for the CSI scheme. FIG. 6(a) shows a 169 pixels in zonal filter 1 , FIG. 6(b) shows a 133 pixels in zonal filter 2, FIG. 6(c) shows a 69 pixels in zonal filter 3, and FIG. 6(d) shows a 25 pixels in zonal filter 4. In order to reduce the computational complexity of this zonal filter 1 , it is proposed next to use the zonal filter 2, shown in FIG. 6(b), the zonal filter 3, shown in FIG. 6(c), and the zonal filter 4, shown in FIG. 6(d). These zonal filters 2, 3 and 4 use zonal masks of 133, 69 and 25 pixels or grid points of _?(/, ,., ) , respectively, to compute each Y in (25) and each _{JUl 2} in (26).

Next for τ= 3, the FCSI schemes described above are developed by the use of zonal filter 2, zonal filter 3 and zonal filter 4. Following the same procedure, used in the previous section, the FCSI scheme for τ= 3 is obtained readily. First the derivations of these FCSI schemes for r = 3 are developed by a means similar to that of the original CSI scheme for r = 2 [10]. Only the terms Y _Mι in (25) and C_{| .,} . in (26) for the FCSI schemes are slightly different from the original CSI scheme. Also in the FCSI scheme, the complexities of the terms Y and A_{} k k} for zonal filter 2, zonal filter 3 and zonal filter 4 are summarized in Tables I and II and are described in the following: 1. By FIG. 6(a), the terms Y_{h h} and A_{hJ k k} for zonal filter 1 are obtained completely in a manner similar to (25) and (26). It uses 169 pixels of R(t t₂ ) to compute each Y ₂ and each A_hJ^_kι . This algorithm is very complex and is just the original CSI scheme.

2. By FIG. 6(b), the terms A_hJιλ^ for zonal filter 2 are obtained, and the result is similar to (26), except that σ = 0 . Each Y in (25) is computed for zonal filter 2 with 133 pixels of ?(., ,., ) . This algorithm still involves extensive computations.

3. By FIG. 6(c), the terms A _kikι for zonal filter 3 are obtained, and the result is similar to (26), but for this case μ = 0 and σ = 0 . Each Y in (25) is computed for zonal filter 3 with 69 pixels of R(t, ,t₂) . This algorithm is also complex.

4. By FIG. 6(d), the terms ,_l7. _,,_.,*, for zonal filter 4 are obtained, and the result is similar to (26), but in this case η = λ = p = μ = σ = Q . Each Y in (25) is computed for zonal filter 4 by the use of only 25 pixels of /?(/, ,._ ) . This algorithm is more compact than any of the other computations given above in cases 1, 2 and 3. TABLE I COMPLEXITY OF A , , _{k k} USING FOUR ZONAL FILTERS.

TABLE II

COMPLEXITY OF COMPUTING Y WITH A COMPRESSION RATIO OF 9 TO 1 FROM AN

IMAGE OF SIZE 512 BY 512 USING FOUR ZONAL FILTERS.

It is shown by computer runs in Section VIII that the zonal filter 4 obtains a PSNR that is similar to any of the other three zonal filters. Thus, this zonal filter 4 represents the most practical and simple zonal filter for the FCSI scheme. The primary advantage of the FCSI scheme with zonal filter 4 over the original CSI scheme is that it substantially reduces the computational complexity. V. CALCULATION OF CONSTANTS

In accord with (27), (28) and the previous section the terms A _kιk2 needed for zonal filter 4 are given by the general formula: if , -A, ≡Omod/?, and 7, -A, ≡Omod/?, if 7, -A, ≡ ±1 mod/?, and ₂ -A₂ ≡ 0mod?₂ or 7, -A, ≡Omod/?, and ₂ -A, ≡±lmod/?, if 7, -A, = ±2 mod/?, and7₂ -A₂ ≡0mod/?₂ or 7, - A, =0 mod /?, and 7, - A₂ ≡ ±2 mod /?, if 7, -A_j ≡±3mod/?, and7₂ -A₂ ≡Omod.?, or 7, -A, ≡ Omod/?, and7₂ -A, = ±3 mod/?,

otherwise or by the array,

The constants , β, γ and δ are the autocorrelation coefficients between the 2-D spline function ?(//?, .??,) In the following analysis, τ and m are assumed to be integers and R(t) is the 1-D cubic-spline function, defined in (1) By the use of (33), one obtains constants α, β, γ and δ as follows

1) Calculation of α By FIG 7 for τ = 3, the value of c. is the sum of the overlapping values of the 2-D spline function .(//?,,/»,) , given by

2) Calculation of β By FIG 8 for τ =3, the value of β is the sum of the overlapping values of the 2-D spline function R{m_x,m₂) and the shifted 2-D splme function R(m_x +r, /.?,), given by

/? = ∑ R(//; /_ ((.--/./)/ r)+ 2 x ∑R((r + //?)/ r)R(/z./r) (35) tιι=\ m=

3) Calculation of γ By FIG 9 for τ = 3, the value of γ is the sum of the overlapping values of the 2-D spline function R(m_x,m₂) and the shifted 2-D splme function -.(/«, +2r,/n.) , given by

7= ∑R{m/τ)R((2τ-m)/τ) (36)

4) Calculation of δ By FIG 10 for τ = 3, the value of δ is the sum of the overlapping values of the 2-D spline function R(m,m₂) and the shifted 2-D splme function /?(/«, + 3τ,m₂) , given by δ = ∑ R((r + //?)/ τ)R((2τ - w)/ r) . (37) m=0 Specifically for τ = 3 the above equations yield the following parameters for α, β, γ and δ: a = l + 4x ∑R²( n/3) = 3.8916 , m=l ? = ∑ R{m /3)«((3 - /«)/ 3)+ 2 x ∑ R((3 + / .)/3)R(//. /3) = 0.3786 ,

/«=! =0 y = ∑R(//?/3 ?((2 x 3 - //i)/3) = -0.1070 , δ = ∑R((3 + /«)/3 ?((2 x 3 - n)/3) = 0.0055 . '"=0 ra=0

FIG. 7 shows the overlap map grid points for calculating α for τ=3, FIG. 8 shows the overlap map grid points for calculating β for τ=3, FIG. 9 shows the overlap map grid points for calculating γ for τ=3, and FIG. 10 shows the overlap map grid points for calculating δ for τ=3.

VI. A NOVEL FCSI ALGORITHM IMPLEMENTED BY WINOGRAD DFT AND

OVERLAP-SAVE METHOD

The Winograd DFT algorithm described in Dean P. Kolba and Thomas W. Parks, "A prime factor FFT algorithm using high-speed convolution," IEEE Trans, on Acoustic, Speech, and Signal Processing, vol. ASSP-25, No.4, pp.281-294, Aug. 1977, [16]; S.

Winograd, "On computing the discrete Fourier transform," Mathematics of Computation, vol. 32, No.141, pp.175- 199, Jan. 1978, [17]; and I. S. Reed, T. K. Truong, R. L. Miller, and B. Benjauthrit, "Further results on fast transforms for decoding Reed-Solomon codes over GF(2ⁿ) for /. = 4, 5, 6, 8," in the Deep Space Network Progress Report 42-50, Jet Propulsion Laboratory, Pasadena, CA, pp.132-155, Jan. and Feb. 1979. [18], the contents of which are hereby incoφorated by reference, with a new type of overlap-save method is used in this section to implement the FCSI scheme for τ= 3.

Consider an image of size N^χN = 512^χ512 pixels. If the FCSI scheme for τ= 3 is used to compress this original image, the size of the compressed image is reduced to |^~512 /3^~|χ |^~512 / 3^~] = 171 x 171 pixels, where

denotes the least integer greater than or equal to x. Since 171 is not a power of two, the 2-D FFT can not be used to solve for the X_k[ in (30). To overcome this problem, one possible method is to expand the data Y in (25) and β in (28), where 0 < 7, , _-, < ?, - 1 for . = 1,2 from 171 >< 171 pixels to

256*256 pixels, by appending zeros to the edges of the original reduced image. Then one would take the 2-D FFT of Y _{h h} and β. _ to obtain Ϋ_{m n} and B_{m n} , respectively. This would yield X_m „ = Ϋ_{m n} I B_{m n} . Finally, one would take the inverse 2-D FFT of the X_{m n} to solve for the X_{k[ kι} , where 0 < A, , A, < 255 , the compressed data to be transmitted or stored.

The disadvantage of using such a 2-D FFT method to compute the CSI scheme for τ= 3 is that it requires more computational time due to the increase in the size of the data from

171 x 171 pixels to 256x256 pixels.

Since the size of the compressed data to be transmitted for τ= 3 is not always equal to a power of two, in order to compress such an odd-sized image, a 9-point Winograd DFT

(WDFT) described below is used to implement the FCSI scheme. First since 171 = 9x 19 the ₂ and β_ . , where 0 ≤ j, ,j₂,s, ,s₂ < 170 , can be divided into 19x 19 sub-images, each of size 9^χ9 pixels, shown in FIG. 1 1. One might think that the 2-D Fourier transform of each sub-image could be implemented directly by the use of the 9 9 WDFT algorithm, and that the 2-D Fourier transforms Y_{m n} and B_{m n} of the Y_jχ and β_ _ of size

171 x 171 could be obtained, respectively. In addition, the X_{m n} could be calculated from the Y_{m n} divided by the B_{m n} .

Finally, one might expect that the inverse 2-D Fourier transform of X_{m n} could be obtained by the use of the 9 9 WDFT algorithm. In other words, the reconstructed image of size 171 x 171 might be obtained. However, this inverse transform is found experimentally to have serious artifacts. This is seen in the reconstructed image, shown in FIG. 12. FIG. 12 illustrates a reconstructed image with serious artifacts that used FCSI method implemented by the direct use of the 9 ^χ 9 Winograd DFT for compression. FIG. 1 1 shows a 19 19 sub-images of size 9 9 in an exemplary image of size 171 x 171 To remove the artifacts found in FIG. 12 by the direct use of the 9x9 WDFT compression described above, a novel type of overlap-save sub-image technique is applied to the FCSI method. An illustrative example of the FCSI algorithm is implemented first by the simplified 5 5 WDFT described below. This new type of overlap-save sub-image method is illustrated in this example, the block diagram of which is depicted in FIG. 13. FIG. 13 is a simple illustrative example of the FCSI algorithm implemented by the

5^χ5 Winograd DFT and the overlap-save sub-image method. The block diagram shown in FIG. 13 is separated into two parts as indicated by broken lines. The first part, labeled by "I", is a FCSI encoder that uses the 5^χ5 WDFT algorithm with the overlap-save sub- image technique. The second part, labeled by "II", is a FCSI decoder. In the first part, consider the source image data to be of size 24x24 pixels as shown in FIG. 13(a). The use of (25) with zonal filter 4 for τ= 3 yields the reduced image needed to obtain the coefficients Y_t of size 8^χ8 pixels as depicted in FIG. 13(b). Next these coefficients Y ■ are divided into four overlapping 5^χ5 sub-images as shown in FIG. 13(c). Note that each sub-image of size 5x5 pixels overlaps each adjacent sub-image with a border of width 2. It is shown experimentally that this border can be used to solve the boundary condition between the two neighboring sub-images when using the Winograd DFT algorithm.

As illustrated in FIG. 13, one uses the next three steps to transform the four 1 overlapping 5^χ5 sub-images of Y_Mι to the corresponding four overlapping 5x5 sub- images of X_k ^ . The first step is to take the 5^χ5 WDFT of the four overlapping 5 5 sub- images of Y _h to obtain the four transformed 5 5 corresponding sub-images of Ϋ_{m n} .

5 The second step is to divide these four 5^χ5 sub-images of Y_{m π} by the 5x5 WDFT of B s_t .s₂ in (33) or β,„ „ for τ= 3 to obtain the four transforms of the 5x5 sub-images of X_m„ . The third step is to take the inverse 5^χ5 WDFT for these four 5^χ5 sub-images of X_{m n} to obtain finally the corresponding four overlapping 5x5 sub-images of X_{k Ji} as shown in

FIG. 13(d). 10 Because some of the pixels in the overlapping border of 5^χ5 sub-images of X_{k k} appear in the other adjacent 5x5 sub-images of X_{kι kι} , the duplicated pixels in the four overlapping 5x5 sub-images of X _{kt kι} are deleted or removed. By this means the four overlapping 5^χ5 sub-images of X_{kι ll2} become the four non-overlapping 4^χ4 sub-images of

X_k . To illustrate this, in FIG. 13(d) each 5x5 sub-image has an overlapping border of

^ width 2; FIG. 13(e) illustrates the remaining samples of each sub-image obtained by the use of this overlap-save method. In FIG. 13(d), first consider all four sub-images with overlapping border in the column direction. Because of the end-effects, the last columns of sub-images 1 and 3 are the duplicated columns in the overlapping border to be removed. ^ However, in sub-images 2 and 4 in FIG. 13(d), the first columns of these two sub- images are also duplicated columns in the overlapping border that need to be deleted. Finally, the above overlap-save method in the row direction is accomplished similarly to that of the column direction. A combination of these four non-overlapping 4^χ4 sub- images of X_kι<k yields the entire 8x8 image of X_{k k} shown in FIG. 13(f). These x_{k k} ⁵ image data are the compressed data to be transmitted or stored. In the second part, using the entire 8^χ8 X_{k kι} compressed data, the 24x24 reconstructed data, shown in FIG. 13(g), are obtained by means of the cubic-spline reconstruction function, given in (32).

Thus the steps included in compressing an image are defining a cubic-spline filter; correlating the filter with the signal to obtain a correlated signal; autocorrelating the filter ^ to obtain autocorrelated filter coefficients; computing a transform of the correlated signal and the autocorrelated filter coefficients; dividing the transform of the correlated signal by the transform of the autocorrelated filter coefficients to obtain a transform of a compressed signal; and computing an inverse transform of the transform of the compressed signal to obtain the compressed signal. ⁵ In order to compress an actual image of size 512x512 by such an overlap-save sub- image method, the 9x9 WDFT instead of the 5^χ5 WDFT is used for the FCSI scheme. Then, the FCSI encoding algorithm using the 9 9 Winograd DFT and the overlap-save sub-image method is summarized in the following steps: • Choose τ= 3. The compression ratio is roughly τ² = 9.

• Apply (25) with zonal filter 4 to the original image of size 512 512 to find all of the

171 x 171 coefficients,' Y 1,1 ,2 . • Also apply (33) to find _{^ s} . Then take the 9x9 WDFT of B_{S t} to obtain the 9x9 coefficients, B,„ „ .

• Divide all of the coefficients Y into the appropriate overlapping 9x9 sub-images of Y _h with a border of width 2.

• Take the 9 9 WDFT of all overlapping 9x9 sub-images of Y to obtain the transforms of the corresponding 9x9 sub-images of Ϋ_{m n} .

• Compute X_{m ιl} = Y_{m n} IB_{m n} to obtain each 9x9 sub-image of X_{m n} .

• Take the inverse 9 9 WDFT of all of the 9*9 sub-images of X_{m n} to obtain the 9x9 sub-images of x_{k k2} .

• Remove the duplicate pixels in the overlapping border of two adjacent 9^χ9 sub-images of X_{kι Aι} . By this means the overlapping sub-images of X_{k ιk} become the non- overlapping sub-images of X_{k k} .

• Combine each non-overlapping sub-image of X_{k k} to obtain the entire image of X_{k k} , where 0 < A,,A, < 511 . These image data X_{k k} are the compressed data to be transmitted or stored. In the FCSI decoding algorithm, the compressed data X_{k kι} constitutes the encoding algorithm; the reconstructed data are obtained by means of the cubic-spline reconstruction function in (32). Hence, by means of the above-described overlap-save sub-image technique, both the compressed and the reconstructed image are obtained and found experimentally to have no apparent artifacts, as shown in FIG. 14. FIG. 14 is a reconstructed image with no apparent artifacts that used the FCSI implemented by the 9x9 Winograd DFT and the overlap-save method for compression. VII. A MODIFIED JPEG ENCODER-DECODER

In this section, a modified JPEG encoder-decoder is presented for image data compression. This algorithm applies the CSI or FCSI scheme with a compression ratio of r² to 1 as the pre-processing step of the JPEG encoder for τ= 2, 3 as shown in FIG. 15. As a consequence the cubic-spline reconstruction with a ratio of 1 to τ² is used for the post-processing step of this JPEG decoder for τ= 2, 3 as shown in FIG. 16. For this compression algorithm, an original image in the RGB (Red, Green and Blue) color space is converted into another preliminary image in YUV color space prior to the CSI or FCSI pre-processing.

This YUV image is followed by the CCIR 601 color space with format 4:1 :1. The size of the original RGB image is assumed to be 512x512x3 =786,432 bytes, i.e., each set of 512x512 = 262,144 bytes is for one of the red, green, and blue colors. After color- space conversion one set of 512x512 bytes is used for Y, and two sets of 256x256=65,536 bytes are used for the U and V color components. It follows from Phillip E. Mattison, Practical Digital Video with Programming Examples in C, John Wiley & Sons, Inc., 1994 [19], the contents of which are hereby incoφorated be reference, that the equations for the conversion from RGB to YUV are F = 0.2995 + 0.587G + 0.1 145, U = 0.493(5 - Y) = 0.4635 - 0.147.fi - 0.289G and V= 0.877(5 - Y) = 0.6155 - 0.515G - 0.1005.

There are two processing steps in the encoder phase. The first step is the preprocessing that uses the CSI or FCSI scheme with a compression ratio of τ² to 1 for each of Y, U, and V images. In this procedure, the input image is a Y image of size 512x512 bytes, and the output image is an encoded image of size f512/ r χ (^"512/ r^"l bytes, where denotes the least integer greater than or equal to x. For the U and V images the input image has 256x 256 bytes so that the output image to be encoded is [^~256/ r]χ |^~256/ r] bytes. In other words, there are two cases used in this step. For the first case τ= 2, the CSI scheme implemented directly with the FFT algorithm is used for the original Y, U and V images and the output images are 256^χ256 bytes for Y image and 128^χ128 bytes for U and V images.

For the second case τ= 3, the FCSI scheme implemented by the 9x9 WDFT with overlap-save method is used for the original Y, U and V images and the output images are 171 x 171 bytes for Y image and 85x85 bytes for U and V images. At the end of the CSI or FCSI algorithm the three separate Y, U, and V images are combined into one YUV image. The second step is to use the JPEG DCT-based encoding algorithm [9]. The image after this step is called the compressed image. This compressed image has now a very small number of pixels when compared to the original image. The resulting image still has the standard JPEG format. As a consequence, this compressed image can use the standard JPEG decoder, also save on storage and decrease the transfer time for a communication.

Also in the modified JPEG decoder, there are two processes used which are reversed in some of the encoding steps. The first step is the JPEG DCT-based decoding algorithm [9]. After this step, the image file is separated into three separate Y, U, and V images. The second step is the post-processing step that uses the cubic-spline reconstruction with a ratio of 1 to r² for Y, U, and V image. This step uses only the cubic-spline function to reconstruct the image data. After this inteφolation the size of the Y image is therefore converted from f512/ r]^χ f512/ r] , i.e., 256x256 for τ= 2 or 171 x171 for τ= 3, to 512 x 512 bytes, and the U and V images are increased from |^"256/ r χ 256/ r] , i.e., 128x 128 for τ= 2 or 85*85 for τ= 3, to 256 x 256 bytes. Then the three Y, U, and V images are combined again into one YUV format. Finally, this YUV image is converted into the reconstructed RGB image. Again, it follows from [19] that the equations for the conversion from YUV to RGB are R = Y + 1.140 , G = Y- 0.395(7- 0.581 V, and 5 = Y +

2.032(7

VIII. EXPERIMENTAL RESULTS

Let X{ι,j) and S(ι,j) be the original and reconstructed images, respectively, where 0 < i < M - 1 and 0 < j ≤ N - 1 are the index numbers separated in the vertical and horizontal directions of the images. The mean-square error (MSE) of the 2-D signal of an image is given by

Thus, the PSNR of the 2-D signal are defined by ^f 255^{2 Λ}

5SNR _v (_Vj5) = 101og_l0 for X e { Y, U, V }, (39) MSE _j where MSE_y , MSE_V and MSE_V are the MSE of color component Y, U and V, respectively The total PSΝR of the color component Y, [/ and V is defined by

255'

PSNR_T(dB) = 10 log 10 ^' (40)

(MSE_Y + MSE_υ + MSE_V ) / 3 Experimental results for the 2-D signal are presented using linear inteφolation, linear spline inteφolation, cubic-convolution inteφolation, cubic B-spline inteφolation and the CSI scheme. Using the encoding method descπbed in Section II and the decoding method descπbed in Section III for the 2-D signal with x= 2, from (38) and (39), the PSΝR values of the 2-D signal with a compression ratio of 4: 1 ( τ - 2 ) are calculated for some gray images of size 512 by 512. The experimental results of the 2-D signal using the above five inteφolation schemes are shown in Table III. The results in this table show that the CSI scheme obtains the best PSΝR over all other compared inteφolation methods.

TABLE III THE PSΝR (dB) OF 2-D FUNCTION WITH COMPRESSION RATIO 4 1 ( r = 2 )

In Table IV, the expeπmental results with a compression ratio of 9: 1 ( r = 3 ) are presented for the FCSI algorithm that uses zonal filter 1, zonal filter 2, zonal filter 3 and zonal filter 4 for the same gray images of size 512 by 512. It is observed from this table that the FCSI scheme with zonal filter 4 obtains a PSNR that is similar to the other three zonal filters. In addition, the experimental results with a compression ratio of 9:1 (r = 3) that use linear inteφolation, linear-spline inteφolation, cubic-convolution inteφolation and the FCSI scheme with zonal filter 4 are shown in Table V for the same gray images of size 12 by 512. It follows from this table that the FCSI scheme with zonal filter 4 obtains the best PSNR of the four inteφolation methods.

TABLE IV THE PSNR (dB) WITH COMPRESSION RATIO OF 9: 1 ( τ = 3 ) FOR THE FCSI SCHEME WITH

DIFFERENT ZONAL FILTERS.

TABLE V

THE PSNR (dB) WITH COMPRESSION RATIO OF 9: 1 ( τ = 3 ) FOR FOUR INTERPOLATION

SCHEMES.

Table VI lists the PSNR values of the reconstructed color test (Lena) image of size 512 by 512 at different compression ratios for the JPEG method [13], the AOL algorithm [8], and the modified JPEG encoder-decoder for τ= 2 described in FIG. 15 and FIG. 16, also see [10]. For the same compression ratios the PSNR of the Lena image, obtained by the modified JPEG encoder-decoder for τ= 2 are higher than those of both the JPEG method and the AOL algorithm. Furthermore, using (40), Table VII lists the total PSNR values of the colored reconstructed Lena image of size 512 by 512 at higher compression ratios for the JPEG method, the AOL algorithm, and the modified JPEG encoder-decoder for τ= 2 and τ= 3. For the same compression ratio of 250:1 the total PSNR value of the Lena image, obtained by the modified JPEG encoder-decoder for τ= 3, is higher than the JPEG method, the AOL algorithm and the modified JPEG encoder-decoder for τ= 2.

TABLE VI PSNR (dB) VALUES OF THE RECONSTRUCTED COLOR LENA IMAGE OF SIZE 512 BY 512 AT DIFFERENT COMPRESSION RATIOS FOR JPEG, AOL, AND MODIFIED JPEG ENCODER- DECODER FOR τ= 2.

TABLE VII TOTAL PSNR (dB) VALUES OF THE RECONSTRUCTED COLOR LENA IMAGE OF SIZE 512 BY 512 AT DIFFERENT COMPRESSION RATIOS FOR JPEG, AOL, MODIFIED JPEG ENCODER- DECODER FOR τ= 2 AND τ= 3.

For the gray Lena image of size 512 by 512, the computational times of the CSI and FCSI schemes for τ= 3 were implemented on a 400-MHz Intel Pentium II personal computer using C-code. In the encoder, the FCSI scheme using the 9-point WDFT with overlap-save sub-images requires around 0.15 sec when compared with around 0.57 sec for the CSI using the FFT. Therefore, the FCSI scheme is faster than that of the CSI scheme.

Finally, the JPEG method, the AOL algorithm, and the modified JPEG encoder- decoder for τ= 2 and τ= 3 were also implemented on the same 400-MHz Intel Pentium II personal computer using C-code. The computational time of the color Lena image of size 512 by 512 at the compression ratio of 200:1 for these four algorithms are given in Table VIII. In the encoding and decoding, the modified JPEG encoder-decoder for τ= 3 requires only 0.42 sec and 0.27 sec when compared with 1.13 sec and 0.34 sec for the modified JPEG encoder-decoder for τ= 2, 1.09 sec and 0.30 sec for the AOL algorithm, and 0.80 sec and 0.65 sec for the JPEG method, respectively. It is by this means the computational time of the modified JPEG encoder for τ= 3 requires, respectively, 0.71 sec, 0.67 sec, and 0.38 sec less time than the modified JPEG encoder for τ= 2, the AOL encoder, and the JPEG encoder. Also, the computational time of the modified JPEG decoder for τ= 3 requires, respectively, 0.07 sec, 0.03 sec, and 0.38 sec less time than the modified JPEG decoder for τ= 2, the AOL decoder, and the JPEG decoder.

TABLE VIII THE COMPUTATIONAL TIME (sec) OF THE COLOR LENA IMAGE OF SIZE 512 BY 512 AT THE COMPRESSION RATIO OF 200: 1 IMPLEMENTED ON A 400-MHz INTEL PENTIUM II PERSONAL

COMPUTER.

FIG. 17 shows the reconstructed image of Lena at the same compression ratio of 100: 1 , using the JPEG method, the AOL algorithm, and the modified JPEG encoder- decoder for τ= 2. Clearly the Lena image using the modified JPEG encoder-decoder for τ= 2 indicates a better subjective quality of reconstructed image than both the JPEG method and the AOL algorithm. Also FIG. 18 shows the reconstructed image of Lena at a higher compression ratio of 200: 1 , using the JPEG method, the AOL algorithm, and the modified JPEG encoder-decoder for τ= 2 and τ= 3. In this figure, the Lena image using the modified JPEG encoder-decoder for τ= 3 indicates a better subjective quality of reconstructed images than the JPEG algorithm, the AOL algorithm and the modified JPEG encoder-decoder for τ= 2.

FIG. 17(a) shows an original Lena image, while FIG. 17(b) depicts a reconstructed image by the JPEG method with PSNRy = 30.76 dB, 5SN5f/ = 33.90 dB, and PSNRy = 33.90 dB. FIG. 17(c) Illustrates a reconstructed image by the AOL algorithm with PSNRy = 30.91 dB, PSNR = 35.60 dB, and PSNRy = 35.75 dB; and FIG. 17(d) represents a reconstructed image by the modified JPEG encoder-decoder for τ= 2 , PSNRy = 31.20 dB, PSNRy = 35.66 dB, PSNRy = 35.85 dB.

FIG. 18(a) shows a reconstructed image by the JPEG algorithm with encoding time = 0.80 sec, decoding time = 0.65 sec, and PSNRT = 27.90 dB, while FIG. 18(b) depicts a reconstructed image by the AOL algorithm with encoding time = 1.09 sec, decoding time = 0.30 sec, and PSNRf = 31.06 dB. FIG. 18(c) illustrates a reconstructed image by the modified JPEG encoder-decoder for τ= 2, encoding time =113 sec, decoding time = 0.34 sec, and PSNRp = 3118 dB, and FIG 18(d) represents a reconstructed image by the modified JPEG encoder-decoder for τ= 3, Encoding time = 042 sec, Decoding time = 027sec,5SN5?-=3119 dB

9-POIΝT WINOGRAD DISCRETE FOURIER TRANSFORM

The algorithm developed in Dean P Kolba and Thomas W Parks, "A pπme factor FFT algorithm using high-speed convolution," IEEE Trans on Acoustic, Speech, and Signal Processing, vol ASSP-25, No 4, pp 281-294, Aug. 1977 [16], the contents of which are hereby incoφorated by reference, for computing the 9-point Winograd discrete Fourier transform (DFT) is shown by

X(k)= ∑x{n)W"^k, A=0, ,8, (41) n=0 where W = _e-^(2πl9^ _1S the 9-th root of unity m the complex number field, and j = v^ϊ. Algorithm for the 9-point Winograd DFT α, = x{\) + x(8) , _-₂ = Λ(1) - x(&) , a₃ = x(2) + x(l) , α₄ = x(2) -x{l) , a₅ = x(4) + \(5) , a_β = (4) - x(5) , a_η = x{3) + x(6) , α_g = x(3)-x(6), α₉ =-α, +α₅, α,₀ =α, -α₃, _?,,

α,₅ = α, + α₃ +a₅ ,

"16 ="2 ^~a4 + «6> «17 = *(0) + αi5 + ^al » //?, =019740α₉, m₂ = 056858α₁₀ , WJ₃ = 037111a, , , m = 054253α₁₂ , m₅ = 010026a, ₃,

/?₆ = 044228_7_|4 , m-_j = yΛ₇ , /?g = 086603_7₈, mg

/«,₀ =086603a, ₆ , c, = x(Q) - m_η , c₂ = /?2 - '/?₃ , c₃ = w, + /«3 , c₄ = m_x + m , c₅ = c, +c₂-c₃, c^ = c_x +c₃ +c₄, c =c_x -c₂-c₄, cg =m_A-m₆, c_g =m₅-m₆, c_l0 =m₄-m₅, c,, = c_g +c₉ +/n₈, c,2 =c +c_x0-m%, c,₃ =-c₉ +c₁₀ +/?_g, c,₄ =λ(0) + α₇ -m₉, ^(0) = «|₇, A^r(l)=c₃-7C₁₁, A'(2) = c₆-7C₁₂, _Y(3) = c,₄ -7«ιo_> ^(4) = c_η - jc ,

X(5) = c_η + jc_xi, λ (6) = c.,4 + 7//?₁₀ , (7) = c₆ + 7c₁₂ , A"(8) = c₅ + 7 , , Thus, the 9-pomt Winograd DFT requires only 8 multiplications, 49 additions, and 2 shifts, a substantially smaller number of computations than other known algorithms THE 5-POINT WINOGRAD DISCRETE FOURIER TRANSFORM The algorithm developed in [16] for computing the 5-poιnt Winograd discrete Fouπer transform (DFT) is given by

where W = _e ^{' (2π,5)} is the 5-th root of unity in the complex number field, and ₇ = J-T Algorithm for the 5-poιnt Winograd DFT is given by- α, = Ϊ(1)+Ϊ(4), a₂

α₇ =α, +α₃, α₈ = v(0) + α₇,

//?, = 095106c;₅ , w, = 153884α₂ , /«₃ = 036327α₄ , m₄ = 055902α₆ , m₅ = {α₇ , c, = (0) - m_s , c₂ = c, + / ?₄ , c₃ = c_x - m₄ , c₄ = /?, - m₃ , c₅ = m₂ - //?, ,

A-(O) = α₈ , (l) = c, - 7^'c₄ , X(2) = c₃ - c_s , X(3) = c₃ + yc₅ , (4) = c₂ + 7c₄ .

Thus, the 5-point Winograd DFT requires only 4 multiplications, 17 additions, and 1 shift, a substantially smaller number of computations than other known algorithms.

In this invention, a new CSI scheme based on the least-squares method with the cubic-spline function has been proposed to compress the image data. It is shown that the CSI scheme implemented by the FFT algorithm yields a better PSNR performance than all other inteφolation methods for the reconstructed image. In addition a fast CSI, called FCSI, is developed for image compression. Such a FCSI scheme requires fewer additions and multiplications in the decimation process than the original CSI scheme. In the FCSI, a fast 9-point WDFT algorithm is used to compute the CSI scheme with zonal filter 4 and an overlap-save sub-image technique is utilized to solve the serious boundary artifacts between any two neighboring sub-images of the actual image. It is shown by computer runs that the FCSI scheme requires 0.42 sec less time than the original CSI scheme for τ = 3.

Finally, this FCSI scheme for τ= 3 is used with the JPEG encoder-decoder to speed up the modified JPEG encoder-decoder for τ= 2 in color image coding. A computer simulation shows that the modified JPEG encoder-decoder for τ= 3 obtains a better subjective quality and PSNR of the reconstructed image than the JPEG method for high compression ratios. Also, it requires less computational time than the JPEG method, the AOL algorithm and the modified JPEG encoder-decoder for τ= 2 in both the encoding and decoding processes.

It is understood that the exemplary schemes and the respective implementations described herein and shown in the drawings represents only exemplary embodiments of the present invention. Indeed, various modifications and additions may be made to such embodiment without departing from the spirit and scope of the invention. For example, the present invention can be implemented utilizing a computer program, a dedicated electronic circuit, or firmware for an electronic image processing chip. Furthermore, those skilled in the art will appreciate that the method of the present invention is applicable to images having various different sizes with different values for τ and different sizes for Winograd DFT. Thus, the description herein of a 9 x 9, or 5 x 5 of Winograd DFT is by way of illustration only, and not by way of limitation. Further, those skilled in the art will appreciate that various different VLSI implementations of the present invention are contemplated. Thus, these and other modifications and additions may be obvious to those skilled in the art and may be implemented to adapt the present invention for use in a variety of different applications including streaming of video signals, combining with MPEG IV standard for still pictures, applying to JPEG 2000 standard, reducing the size of digital pictures (such as in digital cameras), and the like.

Claims

WHAT IS CLAIMED IS:

1. A method performed by a computer for coding a 1-D signal comprising the steps of: defining a 1-D cubic-spline filter by

(3/2)|.|³ -(5/2)|/|²+ 1 ,0<|t|<l

Λ( H -(^{1 2})H³⁺(^{5 2}) f ^~4H+2 ,l<|t|<2 (1)

0 ,2<|t| applying the filter to an input signal Xm with y. = ∑X(m ₊ jτ)R(m) (9) m=-2r+l to compute Y ; applying A_lk = ^"∑ R( + {j-k)τ)R{m)

_∑2r- __2r+, {R(m)Y =a

j-k≡ Omod ?

∑^{2r Λ}__2τ+iR{m + τ)R{m) = β if 7 -A ≡ ±1 mod/? (10) γ2r- !l_2τ+xR(m + 2τ)R(m)=γ if 7-A ≡ +2 mod/? ∑_m=^"-_2r+ι ^Λ(« + 3r)Λ(/n) = if 7 - A ≡ ±3 mod n 0 otherwise where _ _λ can be expressed in circulant form by letting

and where (k - y^')_n denotes the residue (A - j) modulo n , and

B₀= ,B_i=β,B₂ =γ,B₃=δ,B₄=0, -, 5„_₄ =0,5„_₃ = δ ,B_n_₂ =γ , (12) 5_n__| = /? , to compute 5. ; computing the FFT of Y_] and 5 to obtain Y_m and 5_m , respectively; computing X_m =f_M/5_m;and computing the inverse FFT of _?„, to obtain a compressed data X_k .

2. The method of claim 1 further comprising the step of: ιι-\ applying X_k in ∑A_jkX_k = . for O≤ ≤n -1 , (5) to perform a convolution of

l

X_k and R(t) given by R(t) = 2 (1)

to obtain S(t) given by S(t) = ∑X_kΨ_k(t) = ∑X_kR(t-kτ) (3).

*=0 k=0

3. The method of claim 1, wherein τ= 2.

4. The method of claim 1, wherein τ= 3.

5. An encoder for data compression comprising: a 1-D cubic-spline filter defined by

(3/2)|.|³ -(5/2)|/|²+ 1 ,0<|t|<l

R(t) = -(l/2)|.| +(5/2)|.|²-4|/| + 2 ,l<|t|<2 (1) 0 ,2<|t| means for applying the filter to an input signal Xm with

2.-1

Y_j= ∑ X(m + fτ)R(m) (9) to compute 7, : means for applying

2r-l

= Σ R(m + (J- ^■k)τ)R(m) m=~2τ+\

∑²,^r=^"-₂r₊ι(Λ(^/«))² =« if 7 - A ≡ 0 mod «

∑ -m²,==--₂2_rr₊+.lΛ ' (/n + r)Λ(/«) = /i if j — k ≡±lmod/? (10) _"2τ-\

= ' Σ;,V r₊ι R('" + 2τ)R(m) =χ if j - k ≡ ±2 mod n Σ,_.,^r=^~-₂. ₊ι R(m + 3r)Λ(/n) = δ if j - k ≡ ±3 mod n 0 otherwise where A_lk can be expressed in circulant form by letting ,, =V.,„ . ⁽ID and where (A - j)_n denotes the residue (k - j) modulo n , and

B₀=a,B_x =β,B₂ =γ,B, =δ,B =0, •••, 5„_₄ =0,5„_₃ = δ ,5„_₂ =γ , (12)

_?„_, = β , to compute B_] ; means for computing the FFT of Y_j and 5 to obtain Ϋ_m and B_m, respectively; means for computing X_m = Y B_m ; and means for computing the inverse FFT of X_m to obtain a compressed data X _k .

6. The encoder of claim 5 further comprising a decoder having means for

»-l applying X_k in ∑A _kX_k -Y_} for 0≤f≤n-l , (5) to perform a convolution of (3/2)||³ -(5/2)|.|²+ 1 ,o<|t|<ι

X_k and R(t) given by R(t) = -(l/2)|.|³+(5/2)|.|²-4|t| + 2 ,l≤|t|<2 (1) 0 ,2<|t|

to obtain S(t) given by S(t) = ∑ ΣX-_kΨ_k(t) = ∑X_kR(t-kτ) . (3) k=Q k=0

1. The encoder of claim 5, wherein τ= 2.

The encoder of claim 5, wherein τ= 3.

9. A method performed by a computer for coding a 2-D signal comprising the steps of: defining a 2-D cubic-spline filter by R(/,,/,) = R(/₁)-R(/,), (17) where R{t_x) and R(/₂) are 1-D cubic-spline functions, respectively; applying the filter to an input signal X(t ,t₂) with periods n_xτ and /.,_-, with

2r-l 2r-l ^γ _h,l2 = Σ Σ X(m_x+jτ,m₂+j₂τ)R(m,m₂), (25) m,=-2r+l m₂=-2.+l to compute Y_} ; pplying A _k, R(m_x+(j -k)τ,m₂+(j₂-k₂)τ)R(m_x,m₂)

y ,2²r^r-^_l , AR(nι,,ιi ))^{2 i}f 7 - A, ≡ 0 mod /?, and 7₂ - A, ≡ 0 mod /?,

±τ,m₂)R(m ,m₂) = β ⁱf s ±lmod/?, and7₂ -A, ≡0mod/?₂ or -A ≡ Omodn, and7₂ -A, ≡ ±lmod/?,

Σ»,W-2r₊l^Λ(W ±2τ,m₂)R(m ,ιn₂) =γ i 7 -A ≡ ±2 mod «, and ]₂-k₂ ≡ 0 mod /?, or 7, -*. ≡ Omodn, and 7, -A, ≡±2mod/?,

Σ,V« -2r₊l ^R('ⁿ ±3τ,m₂)R{m_x,m₂) = δ i 7 -k ≡ +3 mod /?, and 7, - A, ≡ 0 mod /?, or -A ≡ Omod/?, and 7, -A, ≡±3mod/?₂

Σm^~L,-2τ_*} ^R(^m ± r, //?, ± τ)R{m , //?, ) = r; i 7 -A ≡ ±1 mod/?, and7₂ -A, ≡ ±lmod/?, Σ,.VV-2.₊l ^R(m + 2r, ιn₂ ± τ)R(m_x ,//?.) = Λ if 7 -A ≡±2mod/?, and7₂ -A₂ ≡±lmod/?. or 7 -A ≡±lmod/?, ^ancl7 ^_^₂ ≡±2mod/?,

.R(m ± 3r, //?, ± τ)R{m_x ,m₂) = μ if 7 -A ≡ ±3 mod/?, and7₂ -A₂ ≡ ±lmod/?, or 7 -A ≡±lmod/?, and7₂ -A, ≡±3mod/?,

±2r, //?, ±2τ)R{m_x,m₂) = p if -A ≡±2mod/., and7₂ -A, ≡±2mod«,

±3τ,m₂±2τ)R{m ,m₂) = σ if 7 -A s±3mod/?_] and 7, -A, ≡±2mod/ι, or7, -A, ≡±2mod/__| and7₂ -A, ≡±3mod/?, otherwise where, the array A_{JUι k<k} can be expressed in 2-D circulant form as

^A _Jι , ^{= B}(k,-_Jl)„_ι (*,-.,)„, > (27) and where (A, -j,)„ denotes the residue (k, -

1,2 and

where 0 < s, < «, - 1 for . = 1,2 , to compute 5.. ; computing the 2-D FFT of Y_ΛJ2 and B_Sj to obtain Y_mn and B_mn, respectively; computing X_m „ = Ϋ_mn IB_mn ; and computing the inverse FFT of X_mn to obtain a compressed image X_k .

10. The method of claim 9 further comprising the step of:

Applying X_k^ in ∑ ∑ X_k A _χkι =Y , 0≤j, ≤n, -1 and = 1,2. (21)

*,=0 *,=0 to perform a 2-D convolution of X_{k kl} and R(t_x,t₂) given by R(t, , . , ) = R(/, ) • R(t, ) to obtain a 2-D reconstructed image S(t, , t₂ ) .

11. The method of claim 9, wherein σ = 0.

12. The method of claim 9, wherein σ = μ = 0.

13. The method of claim 9, wherein η = λ = p = μ = σ = 0

14. The method of claim 9, wherein τ= 2.

15. The method of claim 9, wherein τ= 3. 8743

16. An encoder for data compression comprising: a 2-D cubic-spline filter defined by R(t,t₂) = R(t_x)- R(t₂), (17) where R(/,) and R(t₂) are 1-D cubic-spline functions, respectively; means for applying the filter to an input signal X(t_x,t₂) with periods n,τ and

2r-l 2r-l

/?, r , with Y_{h h} = ∑ ∑ X(m_x + j_x τ, //?, + 7, τ)R(m ,m₂), (25) to obtain Y, means for applying h R(m,+ (7, - A, )r, //., + (j₂ - k₂ )τ)R(m , m₂ )

if 7, -k ≡Omod/., andy₂ -k₂ =0mod«₂ if y, -/e, =±lmod«, andy, -/e, ≡Omod/., or y, - , ≡Omodn, andy^", -k₂ ≡±lmod/., if 7, -k s+2mod«, and/, -k₂ ≡O mod/2, or 7, -k_x ≡0 mod/7, andy₂ - , ≡±2mod/z, if 7, - , ≡ ±3 mod2, andy₂ - , ≡Omodfl, or , -/, ≡O mod/2, andy, -k₂ ≡ ±3 mod/2, ify, -k_x ≡±l mod/2, andy, -k₂ ≡±l mod/2. ify, - , 2+2mod«, andy^", -/, ≡±lmod/., ory, - , ≡±l mod/2, andy₂ -k₂ ≡ ±2 mod/2, ify, - , ≡ ±3 mod/2, andy^", -k₂ ≡±l mod/2, ory, - , ≡±l mod2, andy, -A, ≡ ±3 mod/2, if y, - , ≡ ±2 mod.., andy₂ - /, ≡ ±2 mod/., if y, - , ≡ ±3 mod/2, andy₂ - ₂ ≡ ±2 mod/., ory, - , ≡ ±2 mod/i, andy₂ - ₂ ≡ ±3 mod/2,

otherwise

where, the array A _{hhJ k}^_k can be expressed in 2-D circulant form as

^■v.). (27) and where (A, -j,)„ denotes the residue (A, - 7, ) modulo «, for . = 1,2 and ^λ

where 0 <s <n - 1 for / = 1,2 , to obtain B, means for computing the 2-D FFT of Y, , and B . to obtain Y_m „ and B„ respectively; means for computing X_mι, = Y_mn /B_mπ ; and means for computing the inverse FFT of X_mn to obtain a compressed image

17. The encoder of claim 16 further comprising a decoder having means for applying X _k^ in ∑ ∑ X_k^₂ A_{J M,} = Y_{Jι 2} , O≤ , <n, -1 and . =1,2. (21)

*,=0 *₂=0 to perform a 2-D convolution of X_{k k} and R(t_x,t₂) given by R(t,t₂) = R(t)-R(t₂) to obtain a 2-D reconstructed image S(t, , t₂ ) .

18. The encoder of claim 16, wherein τ= 2.

19. The encoder of claim 16, wherein τ= 3.

20. The encoder of claim 16, wherein σ = 0.

21. The encoder of claim 16, wherein σ = μ = 0.

22. The encoder of claim 16, wherein η = λ = p- μ = σ = 0.

23. A method performed by a computer for coding an image X(t,t₂) comprising the steps of: applying a zonal filter /?(»?,,.??,) to the image X(t,t₂) where ,m₂+j₂τ)R(m_],m₂), (25)

to compute Y J\ Jl '

2τ-\ 2τ-\ applying ,_ι ._i|A. = ∑ ∑ R(/«, +(7, -k )τ,m₂+(j₂ -k₂)τ)R(m_x,m₂) m,=-2r+l m-,=-2τ+\ i y-.²ni^r,^_ ιι¹ι,=-2 ,τ+\ ,(/?(/», I,'/« ., ) •')•'² f 7 0 mod /?, and 7₂ - A, ≡ 0 mod «,

,2r-l

-2r+l /?(//?, ± r, m₂)R(ιn_l ,//?,) = β ⁱf 7 -A ≡ ±1 mod /?, and 7₂ - A, ≡ 0 mod //, or 7 -A ≡ Omod/?, and 7, -A, ≡±l mod n₂

.r-l

∑ in, in, =-2r+ R(m_x ± 2τ, //?, )R(m_x ,//?,)=/ ⁱf 7 -A ≡±2mod«, and7₂ -A, ≡Omod/?, or 7 -A ≡ Omod/?, and 7, -A₂ ≡±2mod/?,

2r-l

__2-+,Λ(/«ι ±3r, /.?,)/.(/.?,, //?,) = 5 ⁱf 7 -A ≡±3mod/?_| and 7, -A, ≡Omod/?, or 7 -A ≡ Omod/?, and7₂ -A, ≡±3modn₂

-_2τ+\ ^R('ⁿ\ ^±τ<^m2 ±τ)R(m_l,m₂) = η if 7 -A ≡±lmod/?, and7₂ -A₂ ≡±lmod/?,

2r-l^" jm, ιιt₂=-2τ+\ R(m_x ±2r, //?, ±r)Λ(//?,,//? ) = I if 7 -A ≡±2mod/?, andy, -A₂ ≡±l mod/?, or -A ≡ ±1 mod ?, and 7₂ - A, ≡ ±2 mod /?,

-m, .(//?, ±3r, in, ± τ)R(m ,m₂ ) = // ⁱf 7 -A ≡ ±3mod«, and 7, -A₂ ≡±lmod/?, or 7 -k ≡±lmod/?, and7₂ -A₂ ≡±3mod/.,

Σ« *.-2_Γ+I^Λ(^WI ^±2r-^m2 ±2τ)R(m m₂) = p ⁱf 7 -A ≡±2mod«_] and7₂ -A, ≡±2mod/?, ∑ , -₂.₊,^Λ('«. ± r,»«2 ±2r)/^.(m_l, ₂) = σ if 7 -A ≡ ±3 mod/., and 7, -A, ≡ ±2 mod/?, or , -A, ≡ ±2 mod/?, and 7, -A₂ ≡ ±3 mod/?, otherwise where, the array A _kk can be expressed in 2-D circulant form as

^- 2.*l* ~^B ,-M⁾,_h ( -Jlr-„₂ (27) and where (A, -7,),, denotes the residue (A, - 7.) modulo n, for 2 = 1,2 and

where 0 < s, < n, - 1 for /^' = 1,2 , to compute 5. . ; where constants α, β, γ, δ, η, λ, p, μ, and σ are autocorrelation coefficients between zonal filter R(m , m₂ ) ; dividing all of the coefficients _. into respective overlapping sub- images of β with a border of width 2; computing a Winograd discrete Fourier transform (WDFT) of all overlapping sub-images of B to obtain coefficients, B_m „ ; dividing all of the coefficients X into respective overlapping sub-images of Y , , with a border of width 2; computing a WDFT of all overlapping sub-images of Y to obtain transforms of respective sub-images of Ϋ_{m n} ; computing X _m „ = Y_{m n} I B_{m n} to obtain each sub-image of X_{m n} ; computing an inverse WDFT of all of the sub-images of X_{m n} to obtain sub-images of X_k ^ ; removing duplicate pixels in an overlapping border of two adjacent sub- images of X_{k kι} ; and combining each non-overlapping sub-image of X_{k k} to obtain the entire image of X_kΛ .

24. The method of claim 23, wherein τ= 2.

25. The method of claim 23, wherein τ= 3.

26. The method of claim 23, wherein the autocorrelation coefficient between zonal filter σ = 0.

27. The method of claim 23, wherein autocorrelation coefficients between zonal filter σ = μ = 0.

28. The method of claim 23, wherein autocorrelation coefficients between zonal filter η = λ = p = μ = σ = 0.

29. The method of claim 23, wherein the WDFT is a 9 x 9 WDFT given by

X(k)= A=0,---,8, (41) and the sub-images of Y _Mι are 9 x 9.

30. The method of claim 23, wherein the WDFT is a 5 x 5 WDFT given by,

4

X(k)= ∑χ(n)W"^k, A=0,---,4, (42) and the sub-images of Y_{Jι 2} are 5 x 5. ιι=Q

31. The method of claim 23 further comprising the step of: applying X_k^ in ∑ ∑^^ ,^¹ _.,.,,. 0< , ≤/., -1 and i = 1,2. (21)

A,=0 λ,=0 to perform a 2-D convolution of X_k and R(t ,/,) to obtain a 2-D reconstructed image S (/,,.,) .

32. A modified JPEG encoder for image compression comprising: means for applying a zonal filter R(m,m₂) to an image X(t_x,t₂) where

2.-1 2.-1 ^γ _l„_l = Σ Σ X(m_x+jτ,m₂+j₂τ)R(m,m₂), (25) m,=-2τ+\ m₂=-2r+l to compute Y ; means for applying

2r-l 2r-i

^A _ι,_l2._k,_k = Σ Σ R(m_x+(j_x-k)τ,m₂+(j₂-k₂)τ)R(m,m₂)

;.i_|=-2r+l m₃=-2r+l

8743

≡ 0 mod ?, and 7 - A, ≡O mod «₂

≡ ±lmod//, and7₂ -A₂ ≡Omod/?, ≡ 0 mod /?, and 7₂ - A₂ ≡±l mod «, ≡ ±2 mod /?, and 7₂ - A, ≡ 0 mod «₂ ≡ 0 mod /?, and 7₂ - A, ≡ ±2 mod n , ≡±3mod/?, and 7, -A, ≡Omod//, ≡Omod/?, and7₂ -A, ≡±3mod/?₂ ≡+lmodw, and7₂ -A, ≡±lmod/?, 3 ±2 mod/?, and 7, -A, ≡ ±lmod/?, ≡±l mod/?, and7₂ -A₂ ≡±2mod/?₂ ≡±3mod ?, and7₂ -A, ≡±l mod /?₂ ≡ ±lmod/?, and7₂ —k₂ ≡±3modn, ≡±2mod/?, and7₂ -A, ≡±2mod/?,

≡±3mod/?, and 7, -A, ≡±2mod/?₂ or , -A, s±2mod«, and7₂ -A, ≡±3mod/?₂

otherwise where, the array A _kk can be expressed in 2-D circulant form as

J\Jl * kl.k = B ^ykι- ι) _h \^k2 .')„, (27) and where (A, -y,)„ denotes the residue (A, - 7, ) modulo «, for 2=1,2 and

where 0 < s. < n. - 1 for 1 = 1,2 , to compute _3_S| . , where constants a, β, γ, δ, η, λ, p, μ, and σ are autocorrelation coefficients between zonal filter R(m_x,m₂) , means for dividing all of the coefficients into respective overlapping sub-images of B with a border of width 2; means for computing a Winograd discrete Fourier transform (WDFT) of all overlapping sub-images of B to obtain coefficients, B_{m n} ; means for dividing all of the coefficients Y into respective overlapping sub-images of Y, , with a border of width 2; means for computing a WDFT of all overlapping sub-images of Y to obtain transforms of respective sub-images of Y_{m n} ; means for computing X_{m n} = Ϋ_{m n} l B_{m n} to obtain each sub-image of X_{m n} ; means for computing an inverse WDFT of all of the sub-images of X_{m n} to obtain sub-images of X_kχ ; means for removing duplicate pixels in an overlapping border of two adjacent sub-images of x_{k kι} ; and means for combining each non-overlapping sub-image of X_{k k} to obtain the entire image of X_{k k} .

33. The encoder of claim 32, wherein τ= 2.

34. The encoder of claim 32, wherein τ= 3.

35. The encoder of claim 32 wherein the autocorrelation coefficient between zonal filter σ = 0.

36. The encoder of claim 32, wherein autocorrelation coefficients between zonal filter σ = μ = 0.

37. The encoder of claim 32, wherein autocorrelation coefficients between zonal filter η = λ = p = μ - σ = 0 .

38. The encoder of claim 32, wherein the WDFT is a 9 x 9 WDFT given by

8

X(k) = ∑x( ) W^nk , A = 0,--- ,8 , (41) and the sub-images of Y_{J{ j2} are 9 x 9. «=o

39. The encoder of claim 32, wherein the WDFT is a 5 x 5 WDFT given by,

4

X(k) = ∑ Λ-(H) W"^k , A = o, • • • , 4 , (42) and the sub-images of Y_{h h} are 5 x 5.

40. The encoder of claim 32 further comprising a decoder having means for n,-l «, -! applying X _{k Ji2} in ∑ ∑ X_{kχ k}, A_{l JιJl}^ = Y , 0 < j, ≤ n, - \ aaά i = \, 2. to

A, =0 *₂ =0 perform a 2-D convolution of X_{kt ιk} and R(t_x ,t₂ ) to obtain a 2-D reconstructed image.

41. A method performed by a computer for coding a signal comprising the steps of: defining a cubic-spline filter; correlating the filter with the signal to obtain a correlated signal; autocorrelating the filter to obtain autocorrelated filter coefficients; computing a transform of the correlated signal and the autocorrelated filter coefficients; dividing the transform of the correlated signal by the transform of the autocorrelated filter coefficients to obtain a transform of a compressed signal; and computing an inverse transform of the transform of the compressed signal to obtain the compressed signal.

42. The method of claim 41, wherein the signal, the filter, and the transforms are one dimensional.

43. The method of claim 41, wherein the signal, the filter, and the transforms are two dimensional.

44. The method of claim 41 , further comprising the step of computing a convolution of the compressed signal and the filter to obtain a reconstructed signal.

45. The method of claim 41 , wherein the step of computing a transform comprises computing a fast Fourier transform (FFT) and the step of computing an inverse transform comprises computing an inverse FFT.

46. The method of claim 41, wherein the step of computing a transform comprises computing a Winograd discrete Fourier transform (WDFT) with an overlap- save scheme and the step of computing an inverse transform comprises computing an inverse WDFT.

47. The method of claim 41, further comprising the step of defining a zonal filter to simplify the steps of correlating and autocorrelating.

48. A modified JPEG encoder for data compression comprising: a cubic-spline filter; means for correlating the filter with an input signal to obtain a correlated signal; means for autocorrelating the filter to obtain autocorrelated filter coefficients; means for computing a transform of the correlated signal and the autocorrelated filter coefficients; means for dividing the transform of the correlated signal by the transform of the autocorrelated filter coefficients to obtain a transform of a compressed signal; and means for computing an inverse transform of the transform of the compressed signal to obtain the compressed signal.

49. The encoder of claim 48, wherein the signal, the filter, and the transforms are one dimensional.

50. The encoder of claim 48, wherein the signal, the filter, and the transforms are two dimensional.

51. The encoder of claim 48, further comprising a decoder having means for computing a convolution of the compressed signal and the filter to obtain a reconstructed signal.

52. The encoder of claim 48, wherein the means for computing a transform comprises means for computing a fast Fourier transform (FFT) and the means for computing an inverse transform comprises means for computing an inverse FFT.

53. The encoder of claim 48, wherein the means for computing a transform comprises means for computing a Winograd discrete Fourier transform

(WDFT) with an overlap-save scheme and the means for computing an inverse transform comprises means for computing an inverse WDFT.

4. The encoder of claim 48, further comprising means for defining a zonal filter to simplify the means for correlating and the means for autocorrelating.