CROSSREFERENCE TO RELATED APPLICATIONS

[0001]
This application claims the benefit of prior filed U.S. provisional Application No. 60/400,326, filed on Aug. 1, 2002, the complete disclosure of which is incorporated fully herein by reference.
STATEMENT OF GOVERNMENTAL INTEREST

[0002] This invention was made with Government support under Contract No. NAG58688 awarded by the National Aeronautics and Space Administration (NASA). The Government has certain rights in the invention.
BACKGROUND OF THE INVENTION

[0003]
1. Field of the Invention

[0004]
This invention relates to the field of data compression.

[0005]
2. Description of the Related Art

[0006]
With the explosion of the digital age, there has been a tremendous increase in the amount of data being transmitted from one point to another. Data may be traveling within an office, within a country, from one country to another, from Earth to locations in outer space, or from locations in outer space back to Earth.

[0007]
Increasingly capable instruments and ever more ambitious scientific objectives produce evergreater data flows to scientists, and this is particularly true for scientific missions involving spacecraft. A large variety of compression algorithms have been developed for imaging data, yet little attention has been given to the large amount of data acquired by many other types of instruments. In particular, radar sounders, radar synthetic aperture mappers, mass spectrometers, and other such instruments have become increasingly important to space missions. Although the volume of scientific data obtained has grown with the increasing sophistication of the instruments used to obtain the data, spacecraft capabilities for telecommunications bandwidth have not grown at the same rate. The tightening constraints on spacecraft cost, mass, power, and size, limit the amount of resources that can be devoted to relaying the science data from the spacecraft to the ground. Competition for use of ground station resources, such as the NASA Deep Space Network, further limit the number of bits that can be transmitted to Earth.

[0008]
One approach to increase the “scientific return” in the face of these constraints is to use data compression, as has been adopted for many NASA scientific missions. For example, the Galileo mission to explore the planet Jupiter and its moons has made extensive use of lossy image compression methods, such as the discrete cosine transform, after the high gain antenna of the Galileo spacecraft failed to deploy properly. By compressing the data, the Galileo team was able to capture data using the spacecraft's smaller, properly functioning antenna.

[0009]
Other missions, like NEAR, make routine use of both lossless and lossy image compression to reduce data volume, employing several different algorithms. In both the NEAR and Galileo programs, scientists felt that the inevitable loss of information associated with data compression and decompression was more than compensated by the opportunity to return more measurements, that is, there is net scientific gain when more measurements are returned (or higher temporal/spatial/spectral resolution is achieved), even with loss of fidelity of data returned.

[0010]
Standard image compression methods like discrete cosine transforms and related methods are optimized for image data and are not easily adaptable to the data streams from nonimage sources (e.g., a spectrometer) or to time series data sources (those with a time component, such as video), and their performance characteristics (in terms of what information is lost by compression) are not necessarily optimal for such time series data. One reason for this is that image compression methods take advantage of 2dimensional spatial correlations generally present in images, but such correlations are absent or qualitatively different in timeseries data, such as data from a spectrometer or particle/photon counter. However, the need for compression of nonimage data is growing and will continue to grow in the future. For example, hyperspectral images from a scanning spectrograph are particularly high bandwidth but not suited for compression by standard techniques. Further, lossless compression methods such as Huffman encoding, runlength encoding, and Fast and Rico algorithms, and lossy methods such as straight quantization, provide relatively small compression rates. Thus, it would be desirable to significantly increase the time resolution of such an instrument within the bandwidth allocation currently available, and increase the compression ratio available when compressing this data, while still retaining the scientific value of the compressed data, and while being able to use the same compression method for single or multidimensional applications.
SUMMARY OF THE INVENTION

[0011]
The present invention is a method, system, and computer program product for implementation of capable, general purpose compression that can be engaged “on the fly”. This invention is applicable to the compression of any data type, including timeseries data, and has particular practical application on board spacecraft, or similar situations where cost, size and/or power limitations are prevalent, although it is not limited to such applications. It is also particularly applicable to the compression of serial data streams and works in one dimension for timeseries data; in two dimensions for image data; and in three dimensions for image cube data. The original input data is approximated by Chebyshev polynomials, achieving very high compression ratios on serial data streams with minimal loss of scientific information.
BRIEF DESCRIPTION OF THE DRAWINGS

[0012]
[0012]FIG. 1 is a flowchart illustrating the basic steps performed in accordance with the present invention;

[0013]
[0013]FIG. 2 is a block diagram of the steps of FIG. 1, containing a small sample of simulated timeseries data; and

[0014]
[0014]FIG. 3 illustrates the result of applying the Chebyshev transform of a sample data set.
DESCRIPTION OF THE PREFERRED EMBODIMENT(S)

[0015]
An embodiment of the present invention is described herein with reference to compression of data obtained by a spacecraft, however, it is understood that the claims of the present invention are intended to cover data compression of all kinds of data, for example, medical imaging data, data acquisition/logging applications, and the like. From a technical standpoint, a general purpose compression method for use onboard spacecraft should have the following properties: low computational burden, particularly in the compression step that occurs onboard the spacecraft where there is limited power and CPU cycles; flexibility to work in small data blocks without requiring any specific data block size; minimal demands on buffer memory to perform compression of real time data; compression and decompression of a data block completely independent of that in other blocks so there is minimum loss of information in the event of data dropouts or corruption; and high quantitative fidelity of decompressed data to the original data stream as measured by scientific information content. Use of Chebyshev polynomials as described herein meets all of these criteria.

[0016]
[0016]FIG. 1 is a flowchart illustrating the basic steps performed in accordance with the present invention and FIG. 2 is a block diagram of the same steps, containing a small sample of simulated timeseries data. At step 100, the original input data is divided into blocks to form a matrix of a predetermined size. For onedimensional data, the size of each matrix will be N×1 (i.e., they will have a single “depth” dimension), while for twodimensional applications, it is convenient, though not required to have the data matrices be square, i.e., N×N blocks and for threedimensional applications it is convenient to have the data matrix be a code, i.e., N×N×N blocks. The optimal block size (the value of N) is a compromise among several factors and is not necessarily related to the bit depth. The “best” matrix size for a given application depends on the available computational resources and the nature of the data (i.e., what degree and types of information loss can be tolerated). A larger matrix size can give higher compression performance but will require more computation. Applicant has found N=8 and N=16 to be common acceptable choices though other choices are also acceptable.

[0017]
For applications of the method in two or more dimensions, the original dataset to be compressed actually consists of “frames” that are sampled at a sequence of times. Each frame is an array of data values, in one or more dimensions. An example of a onedimensional array making up a data frame would be the data from a linescan imager (also called a whiskbroom imager), where the dimension in the data array is spatial. Another such example would be the data from a multichannel particle analyzer, where the dimension in the data array could be particle energy. Successive frames (a total of N) would be buffered and interleaved to produce N×M blocks, in which the first dimension is temporal and the other is that corresponding to the data frame of M samples. A second class of applications would be to datasets where each frame is a two dimensional array, such as for video imaging (two spatial dimensions) or for spectral imaging (one dimension spatial and one spectral). To apply the twodimensional compression to this second class of datasets, successive frames also need to be buffered and interleaved, so that N×N blocks, for instance, are formed with time as one dimension. The other dimension is chosen by the user and is, for best performance, the dimension in which the data frames have greater redundancy of information. Since the application illustrated in FIG. 2 is a onedimensional application, each matrix size is 16×1. As can be seen on the left side of FIG. 2, 32 data points are illustrated and they have each been divided into blocks (block B1, block B2) of 16 data points each.

[0018]
At step 102, the Chebyshev transform is applied to a first data block (block B1 in FIG. 2). This results in, in this example, 16 Chebyshev coefficients for the data in block B1. The result of applying the Chebyshev transform to block B1 are shown in the matrix illustrated in FIG. 3.

[0019]
At step 104, for each coefficient in the matrix (the matrix currently being processed), thresholding is performed. This is best illustrated in FIG. 2, where the results of the Chebyshev transform for each block are plotted on graphs G1 and G2. In the example of FIG. 2, a threshold of −10 to +10 has been established, as illustrated by the threshold lines T1 and T2 at these two points (on each graph). In accordance with the present invention, the thresholding step 104 involves the retention of coefficients larger than the given threshold (e.g., in this example, above 10 or below −10). Thus, data points 1, 2, 3, and 5 are retained and the coefficients for data points 4 and 616 are discarded.

[0020]
At step 106, the retained coefficients (those retained after the thresholding process) are quantized. The quantizing is performed because the retained values can be large numbers that would require significant amounts of memory to store; if a floatingpoint processor is used, each number could be as large as 64 bits if stored as doubleprecision. By quantizing the retained values they can be reduced in size to as few as 8 or even 6 bits, depending upon how much compression is required and how much loss can be tolerated.

[0021]
In the example of FIG. 2, the quantization is performed by mapping the amplitudes of the retained coefficients as follows.

[0022]
The basic Chebyshev algorithm applies a quantizer with a fixed step size (uniform quantizer) to the coefficients retained after thresholding. For each block, the maximum and minimum coefficient values, c
_{max }and c
_{min}, are stored and used to determine the quantizer step size. The quantized coefficients Q(i) are calculated as follows:
$Q\ue8a0\left(i\right)=\left({2}^{m}1\right)\ue89e\frac{\left(c\ue8a0\left(i\right){c}_{\mathrm{min}}\right)}{\left({c}_{\mathrm{max}}{c}_{\mathrm{min}}\right)}\ue89e\text{\hspace{1em}}\ue89e\mathrm{for}\ue89e\text{\hspace{1em}}\ue89ei=1,\dots ,\text{\hspace{1em}}\ue89eN$

[0023]
where m is the number of bits to which the coefficients are quantized, N is the block size, and c(i) is the i^{th }retained coefficient in the block. The distortion introduced by uniform quantization can be measured by setting the threshold to zero in the Chebyshev algorithm, forcing the algorithm to retain and quantize all coefficients.

[0024]
Basically, for each plot, the largest and smallest coefficient amplitude is kept and is used to linearly map the other amplitudes into the range 0 to 2^{8 }(0 to 255) where 8 is the number of bits being used to store the coefficients. It is noted that the number of bits can be any number, and the larger the number, the more accurate the reconstructed signal will be, but the lower will be the compression ratio.

[0025]
The Chebyshev coefficients tend to quickly approach zero as j (see equation 5 below) increases. The distribution of coefficient values therefore has a higher mass near zero, and approximating the coefficients better in this high probability region will reduce quantization distortion. This can be achieved by expanding the quantization intervals near zero using a compander function, then uniformly quantizing. A compander compresses dynamic range on encode and expands dynamic range on decode. Ideally the Chebyshev coefficient compander will stretch values near zero while compressing values near the coefficient maximum and minimum.

[0026]
Several compander functions will work to expand the quantization intervals near zero, including logarithmic and trigonometric functions, as well as the Mulaw and Alaw companders generally used for reducing dynamic range in audio signals before quantization. The inverse hyperbolic sine function was found by the applicant to perform particularly well in expanding Chebyshev coefficients near the origin and compressing coefficients away from the origin.

[0027]
If the high frequency coefficients cluster near zero, the compander almost entirely smoothes out the high frequency noise. But along with the advantage of reducing high frequency distortion comes the disadvantage of an often poorly represented DC coefficient. This problem is easily solved by storing the DC coefficient separately and applying the compander only to the AC coefficients.

[0028]
Other known techniques for quantization can be performed, for example, floating point quantization. An example using floatingpoint quantization is discussed in more detail below in connection with a twodimensional example.

[0029]
At step 108, a bit control word is created so that the retained data points can be inserted at the appropriate location in the data signal when it is decompressed, and so that place holders (e.g., zeroes) can be inserted in the appropriate locations where data points have not been retained. In the onedimensional application illustrated in FIG. 2, since the matrix is a 16×1 matrix, there will be 16 bits in each control word; for an N×N, or N×M matrix (where M is the depth of the matrix), there will be N×N or N×M bits in the control words, respectively. In the example of FIG. 2, the control word for the first data block, block B1, is 1110100000000000. This indicates to the decompression program that the first three bits (the first three “ones”) will contain a data bit to be decompressed or reconstructed, the fourth bit (designated by the “0”), will contain a place holder, the fifth bit (“1”) will contain a reconstructed data point, and the remaining bits (all zeroes) will be given place holders.

[0030]
At step 110, the control word created in step 108 can be encoded using lossless compression techniques (an example of which is given in more detail below). By compressing the control word the compression ratio can be significantly increased without any additional data loss.

[0031]
At step 112, a determination is made as to whether or not there are additional data blocks to be processed. Since the data is processed on the fly, data blocks may be continuously accumulating. If there are additional data blocks to be processed at this time, the process proceeds back to step 102, for processing of the next data block. If not, the process proceeds to step 114.

[0032]
At step 114, the compressed data can be transmitted with its control word so that, upon receipt, decompression can take place by, for example, applying the inverse transform as described below.

[0033]
Thus, as noted above, the Chebyshev algorithm is based on three parameters: first, block size, which is the number of samples used per iteration of the compression method; second, the threshold level, the minimum value of the coefficients to be retained; and third, the number of bits, which is the number of bits to which each coefficient is quantized. By varying the parameters for different types of data, good compression ratios can be achieved with minimal error. Higher threshold values always result in better compression ratios at the expense of reconstructed signal quality. Increasing the number of bits generally decreases compression ratio due to the additional bits stored, but gives better reconstructed signal quality. Increasing the block size has more varied results, but a large block size generally increases the compression ration because not as many block maxima and minima are being stored.

[0034]
Following is an example of a preferred embodiment of the present invention used in connection with twodimensional data (e.g., image data).
EXAMPLE

[0035]
Algorithm Outline

[0036]
Onboard the spacecraft, divide the image into square blocks. To each block:

[0037]
Apply the Chebyshev transform

[0038]
Eliminate the lowamplitude transform coefficients

[0039]
Quantize the retained coefficients

[0040]
Encode and transmit the retained coefficients

[0041]
On the ground, decode and apply block artifact reduction algorithm (optional).

[0042]
Onboard Operations

[0043]
Step 1: Block Processing

[0044]
Divide the image into N×N blocks (N=8 in this example). Each block will then have the form
$F=\left[\begin{array}{cccc}{f}_{11}& {f}_{12}& \cdots & {f}_{1\ue89eN}\\ {f}_{21}& {f}_{22}& \cdots & {f}_{2\ue89eN}\\ \vdots & \vdots & \u22f0& \vdots \\ {f}_{\mathrm{N1}}& {f}_{\mathrm{N2}}& \cdots & {f}_{\mathrm{NN}}\end{array}\right]$

[0045]
where f_{ij }is the ij^{th }pixel in the block.

[0046]
This will require buffering N rows of the image before the first block can be processed.

[0047]
If the image size cannot be divided evenly by N, up to N−1 rows and up to N−1 columns need to be padded with the adjacent pixel value.

[0048]
Step 2: Transform

[0049]
Apply the Chebyshev transform to each block:
$\begin{array}{c}{c}_{\mathrm{ij}}=\frac{4}{{N}^{2}}\ue89e\sum _{k=1}^{N}\ue89e\text{\hspace{1em}}\ue89e\sum _{m=1}^{N}\ue89e\text{\hspace{1em}}\ue89e{f}_{\mathrm{mk}}\ue89e{T}_{\mathrm{im}}\ue89e{T}_{\mathrm{jk}}\\ \mathrm{where}\\ {T}_{\mathrm{ij}}=\mathrm{cos}\ue8a0\left(\frac{\pi \ue8a0\left(i1\right)\ue89e\left(j\frac{1}{2}\right)}{N}\right),\text{\hspace{1em}}\ue89ei,j=1,\dots \ue89e\text{\hspace{1em}},N.\end{array}$

[0050]
c_{ij }is the ij^{th }Chebyshev transform coefficient and T_{ij }is a cosine lookup table stored onboard the spacecraft.

[0051]
The transform can also be written in matrix form,
$c=\frac{4}{{N}^{2}}\ue89e{\mathrm{TFT}}^{t}.$

[0052]
Step 3: Thresholding

[0053]
Set all coefficients, c
_{ij}, whose absolute value is less than a commanded threshold value to zero. The exception is the first coefficient, c
_{11}, which has the highest amplitude and is always retained. The retained coefficients are defined by
${c}_{\mathrm{ij}}^{\prime}=\{\begin{array}{cc}{c}_{\mathrm{ij}},& i=j=1\\ {c}_{\mathrm{ij}},& \uf603{c}_{\mathrm{ij}}\uf604\ge \mathrm{threshold},\\ 0,& \mathrm{otherwise}\end{array}\ue89e\text{\hspace{1em}}\ue89ei,j=1,\dots \ue89e\text{\hspace{1em}},N$

[0054]
This results in a matrix whose elements are primarily zero. Only the nonzero coefficients are quantized, encoded, and transmitted.

[0055]
Step 4: Quantization

[0056]
Quantize the retained (nonzero) coefficients by rounding the mantissa of the binary representation to P bits (sign bit included), and storing only Q positive bits in the binary exponent (because the threshold value will always be greater than 1.0, the need for negative exponents is eliminated, saving one bit per retained coefficient).

[0057]
The specific steps are:

[0058]
1) Convert each coefficient to its binary representation, retaining the sign bit (0 positive, 1 negative)

[0059]
2) Shift the radix point to the leading 1 and increment the exponent (normalization).

[0060]
3) Retain the sign bit and P−1 bits beyond the radix point, rounding if the next bit is a 1. Do not store the leading 1 in the mantissa. It is assumed and will save us 1 bit.

[0061]
4) Retain the exponent as a Qbit value.

[0062]
5) Convert the Q exponent bits and the P mantissa bits (sign bit leading) to a decimal integer.

[0063]
For example, quantize the floatingpoint value 673.4375 to 8 bits using a mantissa of P=4 bits and an exponent of Q=4 bits. Following the steps above,

[0064]
1) The binary representation is 0 10101000001.0111 (the sign bit leads)

[0065]
2) Normalized: 0 1.01010000010111×2^{9 }

[0066]
3) Retain the 4 bits 0 011 in the mantissa (note that 1.0101 was rounded to obtain 1.011, and the leading 1 was not retained).

[0067]
4) Retain the exponent 9 (1001 in Q=4bit binary)

[0068]
5) Convert the exponent bits followed by the sign bit and the mantissa bits (1001 0 011) to an integer: 147

[0069]
The leading 1 in the normalization is assumed and must be replaced when converting back to floating point. The restored floatingpoint value is then +1.011×2^{9}, or 704.

[0070]
For 12bit planetary images compressed using 8×8 blocks, Q=4 bits are stored in the exponent of all thresholded Chebyshev coefficients. The mantissa of coefficient c_{11 }is rounded to P=7 (6 bits plus sign bit), and the mantissa of all other coefficients is rounded to P=4 bits (3 bits plus sign bit).

[0071]
By placing the exponent in the most significant bits of the quantized coefficient and the mantissa in the least significant bits, more efficient lossless encoding (step 8) is achieved. For the lossless encoding step, the quantized values are treated as integers.

[0072]
Step 5: Control Word

[0073]
A ‘control word’ stores the original matrix location of each coefficient in a block and must be transmitted along with the retained coefficients.

[0074]
The control word cw for any given block is defined by mapping each coefficient to either 0 (if it set to zero after thresholding) or to 1:
$\mathrm{cw}=\{\begin{array}{ccc}0,& {c}_{\mathrm{ij}}^{\prime}=0& i,j=1,\dots \ue89e\text{\hspace{1em}},N\\ 1,& \mathrm{otherwise}& \text{\hspace{1em}}\end{array}.$

[0075]
Step 6: Zigzag Scan

[0076]
After thresholding, most of the nonzero coefficients will be clustered in the upper left hand comer of the coefficient matrix. Map both the coefficient matrix and the control word matrix to a vector using a zigzag scan. This enables more efficient runlength encoding of the control word cw and coefficient matrix c′_{ij}.

[0077]
For an 8×8 matrix, the zigzag mapping is defined by
$\left[\begin{array}{cccc}{c}_{11}& {c}_{12}& \cdots & {c}_{18}\\ {c}_{21}& {c}_{22}& \cdots & {c}_{28}\\ \vdots & \vdots & \u22f0& \vdots \\ {c}_{81}& {c}_{82}& \cdots & {c}_{88}\end{array}\right]\Rightarrow \text{}\ue89e\left[{c}_{11}\ue89e{c}_{12}\ue89e{c}_{21}\ue89e{c}_{31}\ue89e{c}_{22}\ue89e{c}_{13}\ue89e{c}_{14}\ue89e{c}_{23}\ue89e{c}_{32}\ue89e{c}_{41}\ue89e{c}_{51}\ue89e{c}_{42}\ue89e{c}_{33}\ue89e{c}_{24}\ue89e{c}_{15}\ue89e{c}_{16}\ue89e{c}_{25}\ue89e{c}_{34}\ue89e{c}_{43}\ue89e{c}_{52}\ue89e{c}_{61}\ue89e{c}_{71}\ue89e{c}_{62}\ue89e{c}_{53}\ue89e{c}_{44}\ue89e{c}_{35}\ue89e{c}_{26}\ue89e{c}_{17}\ue89e{c}_{18}\ue89e{c}_{27}\ue89e{c}_{36}\ue89e{c}_{45}\ue89e{c}_{54}\ue89e{c}_{63}\ue89e{c}_{72}\ue89e{c}_{81}\ue89e{c}_{82}\ue89e{c}_{73}\ue89e{c}_{64}\ue89e{c}_{55}\ue89e{c}_{46}\ue89e{c}_{37}\ue89e{c}_{28}\ue89e{c}_{38}\ue89e{c}_{47}\ue89e{c}_{56}\ue89e{c}_{65}\ue89e{c}_{74}\ue89e{c}_{83}\ue89e{c}_{84}\ue89e{c}_{75}\ue89e{c}_{66}\ue89e{c}_{57}\ue89e{c}_{48}\ue89e{c}_{58}\ue89e{c}_{67}\ue89e{c}_{76}\ue89e{c}_{85}\ue89e{c}_{86}\ue89e{c}_{77}\ue89e{c}_{68}\ue89e{c}_{78}\ue89e{c}_{87}\ue89e{c}_{88}\right].$

[0078]
For example, the thresholded coefficient matrix,

[0079]
is mapped to the vector

[0080]
[c_{11 }c_{12 }c_{21 }c_{31 }c_{22 }0 c_{14 }c_{23 }0 0 0 c_{42 }0 0 c_{15 }0 0 0 . . . 0],

[0081]
and the corresponding control word,

[0082]
is mapped to

[0083]
[1 1 1 1 1 0 1 1 0 0 0 1 0 0 1 0 0 0 . . . 0].

[0084]
Step 7: Control Word Encoding

[0085]
To runlength encode the control word vector, it is truncated after the last 1 and a length specifier telling how many bits are in the truncated control word is transmitted. Two extra bits are saved by eliminating the leading and trailing 1s of the truncated control word vector (the control word vector always begins with 1 because c_{11 }is always retained, and it always ends in 1 by definition of the truncation). The length specifier contains log_{2}(N^{2}) bits (6 bits for 8×8 blocks).

[0086]
In the example above, the transmitted coefficients are

[0087]
[c_{11 }c_{12 }c_{21 }c_{31 c} _{22 }c_{14 }c_{23 }c_{42 c} _{15}].

[0088]
The truncated control word vector is

[0089]
[1 1 1 1 1 0 1 1 0 0 0 1 0 0 1].

[0090]
After eliminating the leading and trailing 1s, the control word vector becomes

[0091]
[1 1 1 1 0 1 1 0 0 0 1 0 0].

[0092]
The control word is now 13 bits in length, so the length specifier (in 6bit binary) is

[0093]
[0 0 1 1 0 1].

[0094]
The length specifier and control word are transmitted along with the quantized coefficients.

[0095]
The exception to the abovedescribed encoding scheme is when c_{11 }is the only coefficient retained. In that case, the 64bit control word vector is

[0096]
[1 0 0 0 . . . 0].

[0097]
After truncation and elimination of the leading 1 (which is also the trailing 1), 0 bits are retained in the control word vector, so 0 bits are transmitted. This special case will be encoded with a length specifier of all 1s,

[0098]
[1 1 1 1 1 1].

[0099]
When the decoder sees a length specifier containing all 1s, it will not read a control word and will know that exactly one coefficient was retained in that block, namely c_{11}.

[0100]
Step 8: Coefficient Encoding (Lossless)

[0101]
The quantized coefficients (step 4) are losslessly encoded using DPCM followed by Rice encoding. This is achieved by first grouping together M blocks of the image (M to be determined by transmission packet size and desired compression ratio—the larger the compression ratio and packet size, the larger M). The M blocks of data are processed as in steps 17, quantized coefficients are grouped together based on their block index in preparation for loseless compression. The order of the indices is chosen according to the zigzag scan described above.

[0102]
The quantized coefficients c
^{m} _{ij }(1≦i, j≦8, 1≦m≦M), in the M blocks
$\left[\begin{array}{cccc}{c}_{11}^{1}& {c}_{12}^{1}& \cdots & {c}_{18}^{1}\\ {c}_{21}^{1}& {c}_{22}^{1}& \cdots & {c}_{28}^{1}\\ \vdots & \vdots & \u22f0& \vdots \\ {c}_{81}^{1}& {c}_{82}^{1}& \cdots & {c}_{88}^{1}\end{array}\right],\left[\begin{array}{cccc}{c}_{11}^{2}& {c}_{12}^{2}& \cdots & {c}_{18}^{2}\\ {c}_{21}^{2}& {c}_{22}^{2}& \cdots & {c}_{28}^{2}\\ \vdots & \vdots & \u22f0& \vdots \\ {c}_{81}^{2}& {c}_{82}^{2}& \cdots & {c}_{88}^{2}\end{array}\right],\cdots \ue89e\text{\hspace{1em}},\left[\begin{array}{cccc}{c}_{11}^{M}& {c}_{12}^{M}& \cdots & {c}_{18}^{M}\\ {c}_{21}^{M}& {c}_{22}^{M}& \cdots & {c}_{28}^{M}\\ \vdots & \vdots & \u22f0& \vdots \\ {c}_{81}^{M}& {c}_{82}^{M}& \cdots & {c}_{88}^{M}\end{array}\right]$

[0103]
are grouped as follows:

[0104]
The M 11bit coefficients:

[0105]
[c^{1} _{11}c^{2} _{11 }. . . c^{M} _{11}]

[0106]
The variable number of 8bit coefficients:

[0107]
[c^{1} _{12}c^{2} _{12 }. . . c^{M} _{12 }c^{1} _{21}c^{2} _{21 }. . . c^{M} _{21}c^{1} _{31}c^{2} _{31}c^{2} _{31 }. . . c^{M} _{31 }. . . c^{1} _{88}c^{2} _{88 }. . . c^{M} _{88}].

[0108]
(Note that many of the c^{m} _{ij }are set to zero by thresholding and are not transmitted.)

[0109]
DPCM encoding followed by Rice compression is then applied first to the stream of 11bit coeficients, and then to the stream of 8bit coefficients. This coefficient ordering ensures that the difference between successive coefficients is small enough to make DPCM followed by Rice encoding more efficient than had the coefficients been grouped in some other order.

[0110]
Step 9: Transmission

[0111]
Each group of M blocks is transmitted in one packet. Transmitted first is a packet header. The header might contain information about the location in the image of the starting block, and the number M of blocks in the packet. Next to be transmitted are the M 6bit control word length specifiers, then the M truncated control words, then the losslessly encoded coefficients.
EXAMPLE

[0112]
Transmission of M=3 8×8 blocks:

[0113]
Block 1:

[0114]
After thresholding all coefficients except for c_{11}=10884.0 and c_{21}=−111.085(are zeroed.

[0115]
The quantized coefficients are then

[0116]
[c_{11 }c_{21}]=[10880 −112],

[0117]
or, after converting to binary according to step 4,

[0118]
[1101 0 010101 0110 1 110],

[0119]
which when converted to integers are

[0120]
[1685 110].

[0121]
The control word vector is

[0122]
[1 0 1 0 0 0 . . . 0].

[0123]
After truncation and removal of leading and trailing 1s, the control word vector becomes

[0124]
[0],

[0125]
and the corresponding length specifier is 1, or

[0126]
[0 0 0 0 0 1].

[0127]
Block 2:

[0128]
After thresholding all coefficients except for c_{11}=12019.67, c_{12}=−347.118, and c_{21}=89.045 are zeroed.

[0129]
The quantized coefficients are then

[0130]
[c_{11 }c_{12 }c_{21}]=[12032 −352 88],

[0131]
or, after converting to binary according to step 4,

[0132]
[1101 0 011110 1000 1 011 0110 0 011],

[0133]
which when converted to an integers are

[0134]
[1694 139 99].

[0135]
The control word vector is

[0136]
[1 1 1 0 0 0 . . . 0].

[0137]
After truncation and removal of leading and trailing 1s, the control word vector becomes

[0138]
[1],

[0139]
and the corresponding length specifier is 1, or

[0140]
[0 0 0 0 0 1].

[0141]
Block 3:

[0142]
After thresholding all coefficients except for c_{11}=7932.9 are zeroed.

[0143]
The quantized coefficients are then

[0144]
[c_{11]}=[7936],

[0145]
or, after converting to binary according to step 4,

[0146]
[1100 0 111100],

[0147]
which when converted to an integer is

[0148]
[1596].

[0149]
The control word vector is

[0150]
[1 0 0 0 0 0 . . . 0].

[0151]
This is the special case where no bits are transmitted for the control word and the corresponding length specifier contains all 1s, or

[0152]
[1 1 1 1 1 1].

[0153]
DPCM and Rice encoding is then applied first to the c_{11 }coefficients (which are 11bit values):

[0154]
[1685 1694 1596],

[0155]
then to the remaining coefficients (which are all 8bit values) in the order described in step 8:

[0156]
[139 110 99].

[0157]
This lossless compression results in some binary stream of values:

[0158]
[(c_{11 }binary stream) (remaining coefficient binary stream)].

[0159]
Without the header, the transmission stream for this packet containing M=3 blocks is then
$\underset{\underset{6\ue89e\text{}\ue89e\mathrm{bit}\ue89e\text{\hspace{1em}}\ue89e\mathrm{length}\ue89e\text{\hspace{1em}}\ue89e\mathrm{specifiers}}{\uf613}}{\underset{\underset{\mathrm{Block}\ue89e\text{\hspace{1em}}\ue89e1}{\uf613}}{0\ue89e\text{\hspace{1em}}\ue89e0\ue89e\text{\hspace{1em}}\ue89e0\ue89e\text{\hspace{1em}}\ue89e0\ue89e\text{\hspace{1em}}\ue89e0\ue89e\text{\hspace{1em}}\ue89e1}\ue89e\text{\hspace{1em}}\ue89e\underset{\underset{\mathrm{Block}\ue89e\text{\hspace{1em}}\ue89e2}{\uf613}}{0\ue89e\text{\hspace{1em}}\ue89e0\ue89e\text{\hspace{1em}}\ue89e0\ue89e\text{\hspace{1em}}\ue89e0\ue89e\text{\hspace{1em}}\ue89e0\ue89e\text{\hspace{1em}}\ue89e1}\ue89e\text{\hspace{1em}}\ue89e\underset{\underset{\mathrm{Block}\ue89e\text{\hspace{1em}}\ue89e3}{\uf613}}{1\ue89e\text{\hspace{1em}}\ue89e1\ue89e\text{\hspace{1em}}\ue89e1\ue89e\text{\hspace{1em}}\ue89e1\ue89e\text{\hspace{1em}}\ue89e1\ue89e\text{\hspace{1em}}\ue89e1}}\ue89e\text{\hspace{1em}}\ue89e\underset{\underset{\mathrm{control}\ue89e\text{\hspace{1em}}\ue89e\mathrm{word}}{\uf613}}{\underset{\underset{\mathrm{Block}\ue89e\text{\hspace{1em}}\ue89e1}{\uf613}}{0}\ue89e\text{\hspace{1em}}\ue89e\underset{\underset{\mathrm{Block}\ue89e\text{\hspace{1em}}\ue89e2}{\uf613}}{1}\ue89e\text{\hspace{1em}}\ue89e\underset{\underset{\mathrm{Block}\ue89e\text{\hspace{1em}}\ue89e3}{\uf613}}{}}\ue89e\text{\hspace{1em}}\ue89e\underset{\underset{\underset{\underset{\underset{\mathrm{coefficients}}{\mathrm{transform}}}{\mathrm{compressed}}}{\mathrm{losslessly}}}{\uf613}}{\left[\mathrm{binary}\ue89e\text{\hspace{1em}}\ue89e\mathrm{stream}\right]}$

[0160]
Ground Operations

[0161]
Decoding

[0162]
Decoding occurs on a blockbyblock basis. The coefficients are read from the transmission stream and decoded. The decompression is achieved by applying the inverse Chebyshev transform to the reconstructed coefficient matrix,
$\begin{array}{cc}\text{\hspace{1em}}& {g}_{\mathrm{ij}}=\frac{1}{4}\ue89e{c}_{11}\frac{1}{2}\ue89e\left(\sum _{k=1}^{N}\ue89e{c}_{\mathrm{k1}}\ue89e{T}_{\mathrm{kj}}+\sum _{m=1}^{N}\ue89e{c}_{1\ue89em}\ue89e{T}_{m\ue89e\text{\hspace{1em}}\ue89ei}\right)+\sum _{m=1}^{N}\ue89e\sum _{k=1}^{N}\ue89e{c}_{k\ue89e\text{\hspace{1em}}\ue89em}\ue89e{T}_{\mathrm{kj}}\ue89e{T}_{m\ue89e\text{\hspace{1em}}\ue89ei}\\ \mathrm{where}& \text{\hspace{1em}}\\ \text{\hspace{1em}}& {T}_{\mathrm{ij}}=\mathrm{cos}\left(\frac{\pi \ue8a0\left(i1\right)\ue89e\left(j\frac{1}{2}\right)}{N}\right),i,j=1,\dots \ue89e\text{\hspace{1em}},N,\end{array}$

[0163]
and g_{ij }are the decompressed pixels.

[0164]
Block Artifact Reduction (Optional)

[0165]
Block artifacts are jumps in the pixel values from one block edge to the adjacent block edge, and are more apparent in highly compressed images. They can be reduced without losing any image detail by adding a gradient correction matrix P_{B} _{ k }to each block B_{k}, where k is the number of blocks in the image. The correction matrix adjusts the pixel values such that the adjacent edges of any two blocks have the same means.

[0166]
The correction matrix for a block is defined by

P _{B} _{ k } =a+bX+cY+dX*X+eY*Y, k=1, . . . , number of blocks

[0167]
where
$X=\left[\begin{array}{cccc}0& 1& \cdots & 7\\ 0& 1& \cdots & 7\\ \vdots & \vdots & \u22f0& \vdots \\ 0& 1& \cdots & 7\end{array}\right],Y=\left[\begin{array}{cccc}0& 0& \cdots & 0\\ 1& 1& \cdots & 1\\ \vdots & \vdots & \u22f0& \vdots \\ 7& 7& \cdots & 7\end{array}\right],$

[0168]
and * represents the elementbyelement product of two matrices.

[0169]
The correction coefficients {a,b,c,d,e} can be found by solving a system of equations formed by applying 5 conditions on the block:

[0170]
1) The average value of P_{B} _{ k }over the block B_{k }is zero:

0={overscore (P)} _{B} _{ k } =a+b{overscore (X)}+c{overscore (Y)}+d{overscore (X*X)}+e{overscore (Y*Y)}=a+3.5b+3.5c+17.5d+17.5e

[0171]
2)5) The mean value of each edge of the block matches the mean of the adjacent block's edge.

[0172]
If we label the edges of a given block B as follows

[0173]
with E_{1 }representing the top row of pixels in the block, E_{2 }the rightmost column of pixels, E_{3 }the bottom row, and E_{4 }the leftmost column, then for each block B we can write 2)5) as

δE _{i} =a+b{overscore (X_{E} _{ i })}+ c{overscore (Y_{E} _{ i })}+ d{overscore ((X*X)_{E} _{ i })}+e{overscore ((Y*Y)_{E} _{ i })},

[0174]
where i=1,2,3,4 and δE_{i }is the average of the difference of the means of adjacent block edges.

[0175]
For example, if we have the following configuration of blocks

[0176]
then, for block B
_{3},
$\begin{array}{c}\delta \ue89e\text{\hspace{1em}}\ue89e{E}_{1}=\frac{1}{2}\ue8a0\left[\frac{1}{N}\ue89e\left(\sum _{j=1}^{N}\ue89e{g}_{1\ue89ej}^{{B}_{3}}\sum _{j=1}^{N}\ue89e{g}_{\mathrm{Nj}}^{{B}_{1}}\right)\right],\\ \delta \ue89e\text{\hspace{1em}}\ue89e{E}_{2}=\frac{1}{2}\ue8a0\left[\frac{1}{N}\ue89e\left(\sum _{j=1}^{N}\ue89e{g}_{\mathrm{jN}}^{{B}_{3}}\sum _{j=1}^{N}\ue89e{g}_{j\ue89e\text{\hspace{1em}}\ue89e1}^{{B}_{4}}\right)\right],\\ \delta \ue89e\text{\hspace{1em}}\ue89e{E}_{3}=\frac{1}{2}\ue8a0\left[\frac{1}{N}\ue89e\left(\sum _{j=1}^{N}\ue89e{g}_{\mathrm{Nj}}^{{B}_{3}}\sum _{j=1}^{N}\ue89e{g}_{1\ue89ej}^{{B}_{5}}\right)\right],\\ \delta \ue89e\text{\hspace{1em}}\ue89e{E}_{4}=\frac{1}{2}\ue8a0\left[\frac{1}{N}\ue89e\left(\sum _{j=1}^{N}\ue89e{g}_{j\ue89e\text{\hspace{1em}}\ue89e1}^{{B}_{3}}\sum _{j=1}^{N}\ue89e{g}_{\mathrm{jN}}^{{B}_{2}}\right)\right],\end{array}$

[0177]
where g_{ij} ^{B} ^{ k }is the ij^{th }reconstructed pixel in block B_{k}.

[0178]
Equations 1)5) can be written as the matrix equation
$\left[\begin{array}{ccccc}1& 3.5& 3.5& 17.5& 17.5\\ 1& 3.5& 0& 17.5& 0\\ 1& 7& 3.5& 49& 17.5\\ 1& 3.5& 7& 17.5& 49\\ 1& 0& 3.5& 0& 17.5\end{array}\right]\ue8a0\left[\begin{array}{c}a\\ b\\ c\\ d\\ e\end{array}\right]=\left[\begin{array}{c}0\\ \delta \ue89e\text{\hspace{1em}}\ue89e{E}_{1}\\ \delta \ue89e\text{\hspace{1em}}\ue89e{E}_{2}\\ \delta \ue89e\text{\hspace{1em}}\ue89e{E}_{3}\\ \delta \ue89e\text{\hspace{1em}}\ue89e{E}_{4}\end{array}\right],$

[0179]
and can be solved for {a,b,c,d,e} by taking the inverse of the matrix. Using {a,b,c,d,e}, the correction matrix P_{B} _{ k }can be calculated. The gradientcorrected pixel value g′_{ij }for block B_{k }is then

g _{iJ} =g′ _{ij} +P _{B} _{ k }where i, j=1, . . . , N.

[0180]
If an edge of a given block lies on the border of the image, δE_{i }for that edge is set to zero.

[0181]
Chebyshev approximation is wellknown; however, for the purpose of simplifying the explanation of the abovedescribed inventive process, it is described in more detail below.

[0182]
The Chebyshev polynomials of the first kind are defined by T
_{n}(x)=cos(n arccos x). These are orthogonal polynomials of degree n on the interval −1≦x≦1, with the weight 1/{square root}{square root over (1−x
^{2})}. They satisfy the continuous orthogonality relation
$\begin{array}{cc}{\int}_{1}^{1}\ue89e\frac{{T}_{i}\ue8a0\left(x\right)\ue89e{T}_{j}\ue8a0\left(x\right)}{\sqrt{1{x}^{2}}}\ue89e\uf74cx=\{\begin{array}{cc}0,& i\ne j\\ \frac{\pi}{2},& i=j\ne 0\\ \pi ,& i=j=0.\end{array}& \mathrm{Equation}\ue89e\text{\hspace{1em}}\ue89e\left(1\right)\end{array}$

[0183]
The polynomial T
_{n}(x) has n zeros on the interval [−1, 1], at
$\begin{array}{cc}{x}_{k}=\mathrm{cos}\left(\frac{\pi \ue8a0\left(k\frac{1}{2}\right)}{n}\right)& \mathrm{Equation}\ue89e\text{\hspace{1em}}\ue89e\left(2\right)\end{array}$

[0184]
for k=1, 2, . . . , n. When T
_{m}(x) is evaluated at its m zeros x
_{k }(k−1, . . . , m) given by (2), the polynomials of degree i, j<m also satisfy the discrete orthogonality relation
$\begin{array}{cc}\sum _{k=1}^{m}\ue89e{T}_{i}\ue8a0\left({x}_{k}\right)\ue89e{T}_{j}\ue8a0\left({x}_{k}\right)=\{\begin{array}{cc}0,& i\ne j\\ m/2,& i=j\ne 0\\ m,& i=j=0.\end{array}& \mathrm{Equation}\ue89e\text{\hspace{1em}}\ue89e\left(3\right)\end{array}$

[0185]
The Chebyshev approximation of order N to a function ƒ(x) is defined by an expansion in terms of Chebyshev polynomials,
$\begin{array}{cc}f\ue8a0\left(x\right)\approx \left[\sum _{k=0}^{N1}\ue89e{c}_{k}\ue89e{T}_{k}\ue8a0\left(x\right)\right]\frac{{c}_{0}}{2},& \mathrm{Equation}\ue89e\text{\hspace{1em}}\ue89e\left(4\right)\end{array}$

[0186]
where the N coefficients are given by
$\begin{array}{cc}{c}_{j}=\frac{2}{N}\ue89e\sum _{k=1}^{N}\ue89ef\ue8a0\left({x}_{k}\right)\ue89e{T}_{j}\ue8a0\left({x}_{k}\right)& \mathrm{Equation}\ue89e\text{\hspace{1em}}\ue89e\left(5\right)\end{array}$

[0187]
for j=0, . . . , N−1, and x_{k }are as before the zeros of T_{N}(x).

[0188]
The Chebyshev approximation is exact on the N zeros of the polynomial T_{N}(x) in [−1, 1]. A suitable transformation from Equation (2) enables the N zeroes to give equal intervals along the xaxis. Furthermore, the Chebyshev approximation has an equal error property, whereby the error of the approximation is distributed almost uniformly over the fitting interval, so it is an excellent approximation to the socalled minmax polynomial which has the least value of the maximum deviation from the true function in the fitting interval. Most important, the Chebyshev approximation can be truncated to a polynomial of much lower degree that retains the equal error property. These mathematical properties of the Chebyshev polynomials are key to their usefulness for approximation and interpolation, as well as for compression of time series data.

[0189]
The Chebyshev compression algorithm is a form of transform encoding applied to data blocks. This category of algorithm takes a block of data (1 or 2dimensional) and performs a unitary transform, quantizes the resulting coefficients, and transmits them. Many types of transforms have been used for image compression, including Fourier series (FFT), Discrete Cosine (DCT), and wavelet transforms.

[0190]
The Chebyshev polynomials are related to the DCT as shown in Equation (2). The Chebyshev approximation, because of the equal error and minmax properties, is particularly suitable for compression of timeseries data, independent of correlations.

[0191]
The actual implementation of the Chebyshev approximation to achieve the present invention is numerically simple and the computational burden is low as shown by Equation (5). It is necessary only to calculate linear combinations of data samples with a small number of known, fixed coefficients. Because the coefficients T_{j}(x_{k}) given in Equation (5) are known, they need not be calculated by the processor, but rather are stored as a table lookup or loaded into memory in advance.

[0192]
The Chebyshev approximation in accordance with the present invention can be programmed using the high level language IDL, although other languages can also be used. By way of example, an IDLlanguage embodiment is described. A compression routine was written in which three parameters could be modified to evaluate the technique. These included the following:

[0193]
Block Size The number of samples of the serial data stream during one iteration of the compression method. The routine continues to process “blocks” of the raw data until the entire data set has been compressed.

[0194]
Threshold An adjustable parameter to balance high compression factors without excessive loss of precision. The threshold value will determine the number of coefficients retained within a block of the serial data stream.

[0195]
Bits The number of bits returned from each retained Chebyshev coefficient.

[0196]
To provide some insight into the Chebyshev method of the present invention, consider a serial data set of N measurements taken in blocks of m samples. Initially, coefficients for all m samples within a block are calculated. Because the Chebyshev approximation is exact on the m zeros of the polynomial, retaining all of these m coefficients preserves the original data within the block exactly (within roundoff errors). By thresholding the coefficients, it is possible to reduce the number of coefficients needed to reconstruct the original m samples with sufficient accuracy to be scientifically useful. However, it is necessary to record which of the m coefficients within a block have been retained (as well as the coefficients themselves) in order to reconstruct the data accurately. It is more efficient to use larger block sizes, but there is a tradeoff with accuracy.

[0197]
For JPEG or similar algorithms, a variance analysis is done on a representative data set and the bit allocation is fixed for a particular type of data. A key advantage of the Chebyshev approximation of the present invention is that simple thresholding of the coefficients can be computed in real time on the actual data. This thresholding technique accomplishes the same purpose of the variance analysis in other algorithms because of the equal error and minmax properties of the Chebyshev polynomials.

[0198]
The Chebyshev compression method of the present invention differs from standard lossycompression techniques such as the Graphics Interchange Format (GIF) or the Joint Photographic Experts Group (JPEG) or the Discrete Cosine Transform (DCT) in that the Chebyshev technique does not rely on 2dimensional correlations in a data set. In fact, one of the advantages to the block encoding performed by this method is that it works extremely well on serial (1dimensional) data. The minmax and equal error properties ensure high fidelity in the reconstructed data with absolute errors roughly independent of the signal strength. This has the consequence that the data compression tends to suppress noise but still reproduces the high signaltonoise peaks in the data with maximum percentage accuracy. Applicants are aware of no comparable compression techniques that are as computationally simple and that work as well on timeseries data. These data sets are typical of instruments used in many areas of space science (for example, hyperspectral data, spectra, particle instruments, magnetometers, etc.).

[0199]
The abovedescribed steps can be implemented using standard wellknown programming techniques. The novelty of the abovedescribed embodiment lies not in the specific programming techniques but in the use of the steps described to achieve the described results. Software programming code which embodies the present invention is typically stored in permanent storage of some type, such as permanent storage of a processor located on board a spacecraft. In a client/server environment, such software programming code may be stored with storage associated with a server. The software programming code may be embodied on any of a variety of known media for use with a data processing system, such as a diskette, or hard drive, or CDROM. The code may be distributed on such media, or may be distributed to users from the memory or storage of one computer system over a network of some type to other computer systems for use by users of such other systems. The techniques and methods for embodying software program code on physical media and/or distributing software code via networks are well known and will not be further discussed herein.

[0200]
It will be understood that each element of the illustrations, and combinations of elements in the illustrations, can be implemented by general and/or special purpose hardwarebased systems that perform the specified functions or steps, or by combinations of general and/or specialpurpose hardware and computer instructions.

[0201]
These program instructions may be provided to a processor to produce a machine, such that the instructions that execute on the processor create means for implementing the functions specified in the illustrations. The computer program instructions may be executed by a processor to cause a series of operational steps to be performed by the processor to produce a computerimplemented process such that the instructions that execute on the processor provide steps for implementing the function specified in the illustrations. Accordingly, the figures herein support combinations of means for performing the specified functions, combinations of steps for performing the specified functions, and program instruction means for performing the specified functions.

[0202]
While there has been described herein the principles of the invention, it is to be understood by those skilled in the art that this description is made only by way of example and not as a limitation to the scope of the invention. For example, while specific implementations in one and twodimensional applications are described in detail herein, threedimensional data can be compressed using the same inventive method and specific implementations for doing so are intended to be covered by the present claims and will be readily apparent to artisans of ordinary skill. Accordingly, it is intended by the appended claims, to cover all modifications of the invention which fall within the true spirit and scope of the invention.