Patent Publication Number: US-7590293-B2

Title: Non-iterative method to restore image from transformed component data

Description:
CLAIM OF PRIORITY 
     This application claims priority under 35 USC 119(e) from U.S. Provisional Application Ser. No. 60/467,862, filed May 5, 2003, which application is incorporated herein by reference. 
    
    
     TECHNICAL FIELD 
     The present invention relates to image data processing, and more particularly concerns restoring an image from a transformed version that includes less than all components of the transformed data. 
     BACKGROUND 
     Digitized images require the storage and communication of large amounts of data. Medical imaging, graphics for games, military reconnaissance, space-based astronomy, and many other applications increasingly strain storage capacity and transmission bandwidth. 
     Although image compression techniques exist, some images, for example computerized tomography (CT) medical images require high resolution, high accuracy, and low contrast. Retrieval of these types of images must be lossless, or nearly so. Even less demanding applications are always open to further data compression. 
     A number of sources produce images in a transform domain, such as the frequency domain. For example, magnetic resonance imaging (MRI) scans, ultrasound scans, and computed tomography/microscopy devices output magnitude/phase data rather than spatial pixels to represent an image. Transform-domain images produce large amounts of data. For example, a spatial image 1000 pixels square, with 8-bit grayscale resolution requires one megabyte of data. A typical double-precision Fourier transform of that image produces about 16 MB of data: 8 MB for the magnitude component and 8 MB for the phase component. The frequency or magnitude component carries some of the information corresponding to the underlying spatial-domain image, and the complementary phase component carries the rest of the information. The ability to restore or reconstruct a spatial-domain image from only one of these transform-domain components could halve the bandwidth and/or time required to transmit the full image in either the spatial or transform domain. Storing only one component could reduce by 50% the space requirements on a disk or other medium. The U.S. Food and Drug Agency (FDA) requires that all medical source data has to be stored. Even when an image is generated initially in the spatial domain, it may be desirable to store or transmit the image in a transform-domain form. 
     Copending commonly assigned U.S. patent application Ser. No. 10/124,547, filed Apr. 17, 2002, demonstrates an iterative technique for restoring a spatial-domain image transmitted or stored as a single component of a transform-domain representation. Unlike previous such techniques, it requires no special conditions in the original image. However, being iterative, it requires processing time and capacity to pursue a number of iterations, making it hundreds or even thousands of times more computation intensive than a non-iterative or closed-form technique. Also, restoration is not lossless; the solution can only approach the pixel values of the original image more or less closely. 
     SUMMARY 
     The present invention offers methods for modifying an original image, communicates a single component of the image in a transform domain, and restores the original image in closed-form, non-iterative processes. The invention encompasses methods, systems, and media. 
    
    
     
       DRAWING 
         FIG. 1  is a block diagram of a system including different aspects of the invention, and includes diagrams of image data processed by the system. 
         FIG. 2  is a flowchart showing an example of a method for modifying an original image for restoration. 
         FIG. 3  shows an example computer program for modifying an original image for restoration. 
         FIG. 4  is a flowchart showing an example of a method for restoring a modified image from a single transform-domain component. 
         FIG. 5  shows a computer program for restoring a modified image from a single transform-domain component. 
         FIGS. 6A-6C  present a demonstration that image restoration is lossless. 
         FIG. 7  is a flowchart of another example method for modifying an original image. 
         FIG. 8  shows another example computer program for modifying an original image. 
         FIG. 9  is a flowchart of another example method for restoring a modified image. 
         FIG. 10  shows another example program for restoring a modified image from a single transform-domain component. 
         FIG. 11  is a flowchart of a further example method for modifying an original image for restoration. 
         FIG. 12  shows a further example program for modifying an original image. 
         FIG. 13  is a flowchart of a further example method for restoring a modified image. 
         FIG. 14  shows a further example program for restoring a modified image. 
         FIG. 15  is a flowchart of an additional example method for modifying one or more original images. 
         FIG. 16  is a flowchart of an additional example method for restoring one or more modified images. 
         FIG. 17  shows an application of the methods of  FIGS. 15 and 16 . 
     
    
    
     DETAILED DESCRIPTION 
       FIG. 1  is a high-level schematic of an embodiment of a system  100  capable of hosting the invention. System  100  may be implemented in hardware, software, or a combination of both. Any of a number of source devices  110 , such as magnetic-resonance imaging (MRI) scanner  111 , ultrasound scanner  112 , or satellite scanner  113  may produce the original image. Satellite water-vapor scanners, for example, produce frequency-domain images. Additionally, the original image may be retrieved from a storage device  114 . A frequency-transform domain original image from one of the input devices  110  is represented as  X (1 . . . K 1 ,1 . . . K 2 ), having both magnitude and phase values 1 through K 1  in a first dimension and 1 through K 2  in a second dimension. Transform image  X  corresponds to an underlying original image in a spatial or other base domain, which is denoted as X(1 . . . N 1 ,1 . . . N 2 ) for an N 1 ×N 2  pixel image. In this example, the base-domain data points must have real values. However, complex-valued data can be treated as two sequences each having real values. 
     Although two-dimensional images are common, system  100  can process images in any number of dimensions. In particular, a two-dimensional image is usually treated as N (N=N 1  or N 2 ) separate one-dimensional pixel sequences X 1 (1 . . . K 1 ) . . . X N (1 . . . K 1 ). Because the Fourier transform is separable, a two-dimensional image can be decomposed by rows or by columns, each row or column processed with a one-dimensional transform, and the N results stored as a two-dimensional transform X(1 . . . K 1 , 1 . . . K 2 ). This method offers considerable computational savings over direct two-dimensional processing. 
     The K 1 ×K 2  values of the transform-domain image components  X  are rational values to produce a lossless restoration. In a digital system, all representable values have a finite precision, and thus are effectively rational. Since all rational values can be considered as integers simply by changing representation scale, the terms “rational” and “integer” are treated as equivalent herein. 
     Processor  120  receives input image  X  over bus  121 , and stores it in a storage facility  122 , which may include internal storage and/or a removable medium for holding image  X  and program instructions for transforming and manipulating the image. Storage  122  also holds a marker quantity  C , which may comprise a single complex value, as explained below. Storage  122  further holds a modified transformed form  XC  of image  X . If the transform is a Fourier transform,  XC  has two complementary components, a magnitude  MXC (1 . . . K 1 ,1 . . . K 2 ) matrix and a phase  PXC (1 . . . K 1 ,1 . . . K 2 ) matrix. Only one of these components, the magnitude  MXC , is output to a communication medium  130 . The complementary phase component  PXC  may be discarded. The term “communication” has a broad meaning herein; device  130  may include a medium for transmitting the image component, and/or for storing it for later restoration. In this embodiment,  MXC  is treated as the modified magnitude component of a discrete Fourier transform (DFT) of a two-dimensional spatial image X(1 . . . N 1 ,1 . . . N 2 ). Each pixel of the spatial image may represent a gray-scale amplitude, one or more color intensities, or any other characteristic, in any desired format. 
     Processor  120  sends the image magnitude component  MXC  to medium  130  for storage and/or transmission to processor  140 . The magnitude component may be further encoded for storage or transmission if desired. An advantage of the present technique is its transparency to any form of lossless compression or encryption. For example, a lossless Lempel-Ziv-Welch (LZW) compression algorithm can be applied to  MXC , so as to achieve an additional 50% reduction of data volume. 
     Processor  140  receives transformed image component after transmission by and/or storage in medium  130 . Processor  140  inverse transforms the communicated component, and removes a marker  C  in the base domain. marker C is equivalent to marker  C  in the transform domain, as explained below. Marker value(s) may be transmitted between the processors along with component  MXC , or may be made known to both processors in other ways. One possible way to communicate the marker value is to embed it within  MXC , in one or more locations where the data points are known to have very small values, such as the high-frequency DFT coefficients. For security, the marker location may be hidden, or it may be scrambled over several locations, making reconstruction difficult without a key. Difficulty of detection is enhanced if the real/imaginary or magnitude/phase values of the marker are embedded so as to maintain symmetry properties of the DFT. 
     Memory  142  stores a base-domain form C of embedded marker  C  along with the final restored base-domain image QF(1 . . . N 1 ,1 . . . N 2 ), as well as program instructions for restoring the image, as in any conventional processor. Processor  140  outputs final image Q on bus  141  to output devices  150 . Examples of output devices include a display  151 , a printer  152 , and a data analyzer for further processing of the restored image. 
     Some applications may desire a transform-domain form of the restored image. Storage  142  includes this form as two arrays representing the magnitude  MQ (1 . . . K 1 ,1 . . . K 2 ) and phase  PQ (1 . . . K 1 ,1 . . . K 2 ) of the DFT-domain restored image. Because the phase component of the original transformed image was discarded before transmission or storage of the (modified) magnitude component, this restoration is known as phase retrieval. 
       FIG. 1  shows elements  110 - 122  as separate from elements  140 - 153 . However, some applications may modify and store an image for restoration on the same computer at a later time. In such cases, only a single processor, storage, and/or bus are employed. Although the marker is normally embedded before communication, it could be embedded after the component is received or retrieved. System  100  may be realized as one or more general-purpose computers, as specialized systems, as elements embedded within other equipment, or in any other physical form appropriate for a particular application. 
       FIG. 2  is a flowchart of an example method  200  for modifying a transform-domain image, representing an image in a spatial or other base domain, for communication as single component in the transform domain. For continuity of exposition, this example assumes a two-dimensional spatial image X(1 . . . N 1 ,1 . . . N 2 ) having pixels expressed as numbers that are rational in terms of a preselected precision. The precision may differ for different embodiments, or for different images in the same embodiment. Standard double-precision floating-point pixel values are assumed for this example, although integer and other floating-point representations are possible. 
     Block  210  receives an image from an image-generating device, a storage, or other input. Normally, this is already expressed in the transform domain, either because it is generated in that form or because it has been converted externally to that form. However, if an application prefers to convert it internally from a base-domain form X(1 . . . N 1 ,1 . . . N 2 ) to transform-domain  X (1 . . . K 1 ,1 . . . K 2 ), block  210  may also perform a conversion by way of a discrete Fourier transform (DFT) or other technique. The number of frequencies K 1  and K 2  in each dimension are selected according to known principles for an appropriate resolution in the restored image. Mathematically,  X (1 . . . K 1 ,1 . . . K 2 )=Re[ X (1 . . . K 1 ,1 . . . K 2 )]+jIm[ X (1 . . . K 1 ,1 . . . K 2 )], the sum of real and imaginary components. These components are also called the Cartesian components of  X . In polar coordinates, the components are called the magnitude  MX (1 . . . K 1 ,1 . . . K 2 ) and the phase  PX (1 . . . K 1 ,1 . . . K 2 ). The real component is complementary to the imaginary component in that together they represent all of the information in the original base-domain image X; the magnitude component complements the phase for the same reason. 
     Block  220  generates a known marker value  C . For most applications, the marker is a single complex number, expressed in polar coordinates as a magnitude (or radius)  MC  and a phase  PC  such that  C = MC exp(j PC ), or equivalently in Cartesian form as the sum of real and imaginary parts,  C = RC +j IC  and is added to every element of  X . That is, the modified transform image  XC (1 . . . K 1 ,1 . . . K 2 )={Re[ X (1 . . . K 1 ,1 . . . K 2 )]+ RC }+j{Im[ X (1 . . . K 1 ,1 . . . K 2 )+ IC ]}, in Cartesian form. It is usually more convenient to set the marker radius as a function of the pixel values in the base-domain image X(1 . . . N 1 ,1 . . . N 2 ). 
     A convenient function for setting the marker radius is  MC =A×N 1 ×N 2 , an amplification factor A times the number of pixels in the underlying base-domain image. The pixel number may represent the size of the underlying base-domain image X, or the size or a largest image to be sent by a system, or in some other manner; these numbers can be approximate or estimated. The amplification factor A may be selected primarily from the maximum pixel magnitude 2^q of the underlying image X. Other possible factors include the precision of the numerical representations for the images, and rounding performed later in the process; the size of A may also be approximate or estimated. A safe lower bound lies in the neighborhood of A=2^q×10^2. Another approximate empirical formula is A=10^5/2^q. Other values are also feasible, up to several orders of magnitude higher, and lower to some extent. The amplification factor may be selected individually for each image, or may be fixed at one value for all images in a set or in a system. Maximum radius is not critical; it need only be low enough so as not to cause ambiguity in the pixel values of the restored image. The marker must be at least large enough to guarantee that the modified image is minimum-phase in the transform domain. Maximum value usually depends upon the available software or hardware floating-point precision. 
     The marker phase  PC  can assume any value except those near 0°, 90°, 180°, or 270°. As a practical matter, using standard floating-point processing, the phase should not lie within the range of ±1° from the four forbidden values. This range can be decreased if desired, by increasing the value of the marker radius  MC .  PC =45° is a convenient approximate value, and may simplify calculation in some instances. 
     For some ordinary applications, the marker  C  remains constant, so that block  220  merely retrieves it from storage for multiple images. Other environments may request block  220  to calculate a separate marker magnitude and/or phase for some or all images. 
     Block  230  embeds marker  C  within image  X  to produce a modified original image  XC . Where  C  has a single complex value, that value is added to every element of the transform image; in the Cartesian formulation,  XC (1 . . . K 1 ,1 . . . K 2 )={Re [ X (1 . . . K 1 ,1 . . . K 2 )]+ RC }+j{Im[ X (1 . . . K 1 ,1 . . . K 2 )]+ IC ]}, where  C  is the sum of its real and imaginary components,  C = RC +j IC . 
     The transform-domain marker in this example represents the transform-domain equivalent of a delta functional in the base domain. Therefore, an base-domain marker C=cδ(n 1 ,n 2 ) could alternatively be applied to a pixel at location X(n 1 ,n 2 ) of the base-domain image, if that image is available at the source. For an image formed as a sequence of discrete pixels, applying a delta functional amounts to adding a predefined constant value to only one of the existing pixels, conveniently to the first element, such as X(1,1) in a two-dimensional image. More generally, it is possible to define markers as functions other than transform-domain constants or base-domain delta functionals. 
     Block  240  calculates the magnitude component of the modified transform-domain image from its real and imaginary parts as  MXC (1 . . . K 1 ,1 . . . K 2  )={[Re XC (1 . . . K 1 ,1 . . . K 2 )^2 +Im XC (1 . . . K 1  . . . 1 . . . K 2 )^2]}^0.5, the square root of sum of the squares. Block  250  sends only the magnitude component  MXC  to medium  130 .  FIG. 1 , for transmission and/or storage. If the marker  C  is redetermined for each image, its values may be sent along with the magnitude as, e.g., part of a file descriptor or header, or within the image data itself, as noted above. 
       FIG. 3  shows a computer program for carrying out a method  300  for modifying a base-domain image for communication as single component in a transform domain. The program is written in the Matlab® language, publicly available from The Math Works, Inc. The original input image may have any number of dimensions. The variable names are explained in the figure. 
       FIG. 4  is a flowchart showing an example of a method  400  for restoring a modified image from a single transform domain component. Block  410  receives the magnitude component  MXC  from medium  130 ,  FIG. 1 . Block  420  performs an inverse DFT or other transform upon  MXC  to produce an approximate restored image QAC(1 . . . N 1 ,1 . . . N 2 ) in the base or spatial domain. 
     Blocks  430  correct for the marker added in methods  200  and  300  above. Block  431  separated QAC into an even part EQAC and an odd part OQAC. Mathematically, EQAC(1 . . . N 1 ,1 . . . N 2 )=[QAC(1 . . . N 1 , 1 . . . 1)]/2. QAC(N 1  . . . 1,N 2  . . . 1) is generated by reversing QAC(1 . . . N 1 ,1 . . . N 2 ) and shifting it circularly one place to the right. The following section gives an example of this procedure; it is a conventional block in most digital signal processor (DSP) integrated circuits. The corresponding odd sequence is OQAC(1 . . . N 1 ,1 . . . N 2 )=[QAC(1 . . . N 1 , 1 . . . 1)]/2. 
     Block  432  subtracts the base domain equivalent form C of the marker from the even part. In this example, C represents a delta function to be subtracted from the first pixel location {1,1}, so block  432  performs EQA(1,1)=EQAC(1,1)− MXC . The even part is then scaled to remove the marker phase “EQA(1,1 . . . N 1 ,1 . . . N 2 )=EQAC(1 . . . N 1 ,1 . . . N 2 )/cos(PC). Because the first element of the odd part is zero when pixels have real values, no subtraction need be undertaken, and the odd part must only be scaled for the marker phase  PC , by OQAC(1 . . . N 1 ,1 . . . N 2 )/sin( PC ). 
     Block  440  resolves the approximate parts to produce exact quantities for restoring the original pixel values from their property of having rational values. The even part is rounded to create an exact or final even part EQ(1 . . . N 1 ,1 . . . N 2 )=Rational[EQA(1 . . . N 1 ,1 . . . N 2 )]. In this example, the function Rational(V)=[Integer(2×10^d×V)]/(2×10^d), where d is the number of decimal places in the representation of its argument V; equivalent functions for binary or other number representations are easily defined. The purpose of this function is to convert a value that lies halfway between two possible values of a pixel in the original image into the correct one of those values; the following section explains why this operation produces exact pixel values. Other functions that achieve this result are also acceptable. Block  440  operates similarly upon the odd part to produce an exact odd part OQ(1 . . . .N 1 ,1 . . . N 2 )=Rational[OQA(1 . . . N 1 ,1 . . . N 2 )]. 
     Block  450  combines the exact even and odd parts to produce the final restored image
 
 Q (1  . . . N   1 ,1  . . . N   2 )= EQ (1  . . . N   1 ,1  . . . N   2 )+ OQ (1  . . . N   1 ,1  . . . N   2 ).
 
Block  460  then outputs Q as the final restored image having pixel values that are the same as those in the original image X(1 . . . N 1 ,1 . . . N 2 ) for a lossless reconstruction.
 
     Some applications may wish to reconstruct the full transform-domain image in addition to—or even instead of—the base-domain image Q. Blocks  470  carry out this task. Block  471  calculates the real part of the restored transform-domain image
 
   RQ   (1  . . . K   1 ,1  . . . K   2 )= DFT[EQ (1  . . . N   1 ,1  . . . N   2 )
 
and the imaginary part
 
   IQ   (1  . . . K   1 ,1  . . . K   2 )= DFT[OQ (1  . . . N   1 ,1  . . . N   2 ).
 
Block  472  computes the retrieved phase component as
 
   PQ   =arctan[   IQ   (1  . . . K   1 ,1  . . . K   2 )/   RQ   (1  . . . K   1 ,1  . . . K   2 ).
 
Because the received magnitude  MXC  contains the embedded marker value, block  472  computes the true magnitude component
 
   MQ =SQRT{ RQ   (1  . . . K   1 ,1  . . . K   2 ) )^2+ IQ (1  . . . K   2 ) )^2}.
 
These calculations could produce polar-coordinate forms instead of the Cartesian forms described. Block  480  outputs the two arrays  RQ  and  IQ  or polar-coordinate magnitude and phase, if desired.
 
       FIG. 5  shows a computer program for carrying out a method  500  for restoring an image communicated as a single component in a transform domain. The program is written in the Matlab® language. The transform image component may have any number of dimensions. The variable names are explained in the figure. The quantity eps in the code is the smallest number representable by the system, to forestall computational artifacts; other implementations may not require it. 
     Explanation of Operation 
     The foregoing describes a representative form of communicating only one component of a transform-domain image, and then restoring the entire image in a closed form—i.e., without iteration. This section demonstrates why the restoration can be lossless. The text refers to numbered equations in  FIG. 6 . Without loss of generality, this section assumes restoration from a magnitude component, one-dimensional sequences, zero-origin indices, and integer data-sequence values, for clarity of exposition. Complex numbers Z are represented interchangeably in the equivalent Cartesian form Z=ReZ+jImZ and polar or magnitude/phase form Z=(MagZ)×exp(−jPhZ), where j=√−1 and “exp” denotes the exponential function. 
     Eqn. 601 defines each term  X (k) of the forward or analysis discrete Fourier transform (DFT) of an N-point sequence having data values X(n), 0≦n≦N−1. Eqn. 602 shows the inverse or synthesis form of the DFT. Each term  X (k) is a complex quantity in the frequency or transform domain, and can be represented in either form of Eqn. 603.  R  and  I  prefixes denote complementary real and imaginary components;  M  and  P  denote complementary magnitude and phase components. 
     A number of transforms, including the DFT, have the linearity property of Eqn. 604, where X 1 , X 2  are any two N-point sequences in a base domain, and a, b are simple scalar constants. 
     A property of any sequence X is that it can be represented as the sum of an even sequence EX and an odd sequence OX, both having N points. Eqn. 605 defines the terms EX(n) of the even sequence. The notation X(−n) signifies the nth term of a sequence that is reversed and circularly shifted one place to the right from the original sequence X. For example, an original sequence {1, 3, 5, 2, 4, 6} is reversed and shifted to {1, 6, 4, 2, 5, 3}, and its even sequence becomes {1, 4.5, 4.5, 2, 4.5, 4.5} according to Eqn. 605. The DFT of any even sequence is entirely real (i.e., not complex-valued), and equals the real part of the full original sequence, as shown in Eqn. 606. Eqn. 607 defines the odd sequence terms OX(n). The odd sequence in the above example is {0, −1.5, 0.5, 0, −0.5, 1.5}. The DFT of any odd sequence is entirely imaginary, and equals the imaginary part of the original sequence, Eqn. 608. 
     Eqns. 605 and 607 demonstrate that, if a sequence has only integer values, then the even and odd sequences consist entirely of integer or half-integer values. Extending this property from integer to rational values reveals that, to any selected scale of a digital number representation, the even and odd sequence values are either exact at that scale or halfway between two exact units. 
     The form of a delta functional or impulse for discrete sequences has an amplitude C for only one of the sequence elements—usually taken to be the first element—and is zero for all other elements. The DFT of a delta in the base domain is a complex constant  C  that is the same for all frequencies k, and has a value given by Eqn. 609. Eqn. 610 shows the polar and Catesian forms of  C .  MC  is the marker magnitude or radius;  PC  is the marker phase. Because of the linearity property of the DFT, Eqn. 604, adding a delta distribution to a base-domain sequence X, Eqn. 611, is is equivalent to adding the constant  C  to every term of the sequence&#39;s DFT  X , Eqn. 612. 
     Eqn. 613 calculates the magnitude component  MXC  for each term of the modified DFT image  XC . This magnitude may be stored or transmitted by itself for lossless restoration of the original base-domain image X, or of the entire transform-domain image  X . 
     Restoration of the original image X from only the component  MXC  employs the DFT properties of even and odd sequences shown in Eqns. 606, 608, and the previously noted fact that a function is the sum of its even and odd parts. Eqn. 614 applies an inverse DFT which, by Eqn. 606, produces a sequence of base-domain terms EQAC(n) that represent an even sequence for restoring the original image. Eqn. 615 produces a sequence of terms OQAC(n) representing a corresponding odd sequence in the base or spatial domain. 
     Sequences EQAC and OQAC still contain a base-domain form C of the transform-domain marker  C  added in Eqn. 612. In this example, the transform-domain marker  C  was applied to correspond to a base-domain delta functional of amplitude C at the first image location; that is, XC(0)=X(0)+C. The marker magnitude only affects the even sequence, and thus need be subtracted only from that part in Eqn. 616. The divisor, cos( PC ), removes the known phase  PC  of the marker. Eqn. 617 removes the marker phase from the odd sequence with divisor sin( PC ). Restricting the marker phase away from the vicinity of the zeros of these two divisors prevents them from introducing overrun or divide-by-zero faults in a digital processor. 
     For sufficiently large marker values, Eqns. 616 and 617 produce even and odd sequences EQA and OQA that are very close to—but not exactly—the even and odd parts of the restored image. The difference from the true sequence-element values arises from the possible introduction of half-integer values into the even and odd sequences by Eqns. 605 and 607. However, the restriction of the original sequence elements to integer values allows the half-integer items to be found and corrected. Eqns. 618 and 619 produce exact, integer even and odd sequences EQ and OQ from the approximate sequences EQA and OQA respectively. The second form of Eqns. 618 and 619 may be used for pixel values having any finite number d of decimal points. 
     Finally, because any function is the sum of its even and odd parts, Eqn. 620 yields the pixel terms of the lossless restored image, Q=X. 
     Restoring the magnitude and phase components then follows from the definitions of those components in Eqn. 603. After transforming the restored base-domain image DFT(Q)= RQ +j IQ , Eqn. 621 calculates terms of the restored magnitude component  MQ = MX , and Eqn. 622 calculates terms of the phase  PQ = PX . 
     Restoration from Phase Component Only 
       FIG. 7  is a flowchart of an example method  700  for modifying an image in a spatial or other base domain for communication as a single phase component in the transform domain, instead of using the magnitude component as in method  200 . Block  711  transforms the image to the transform domain, here again by a discrete Fourier transform. Although the marker could have been added in the base domain, block  720  generates or selects a transform-domain marker of value C having a magnitude R, determined in substantially the same way as previously described in connection with marker  MC  in method  200 ,  FIG. 2 , above. The marker is complex, having a phase arbitrarily set to 45°. Block  730  embeds the marker by adding it to every element of the transformed original image. Block  740  calculates the phase component, and block  750  communicates (transmits or stores) only the phase component as a representation of the entire modified image. 
       FIG. 8  is a Matlab® program  800  for modifying an image for communication as only a phase component in the transform domain. The code sets a marker value C having a magnitude R and a complex phase phi=45°. An array PXIC(1 . . . K 1 ,1 . . . K 2 ), representing the phase elements of the modified image XIC(1 . . . K 1 ,1 . . . K 2 ) may then be communicated—i.e., transmitted and/or stored—on a medium such as 130,  FIG. 1 . 
       FIG. 9  is a flowchart of an example method  800  for recovering or restoring the modified image from phase component  750 ,  FIG. 7 , shown in  FIG. 9  as block  910 . The marker magnitude and phase may also be acquired here if desired. Block  920  forms a complex discrete Fourier transform from the phase component, and block  930  performs an inverse transform, producing a base-domain image. Block  940  removes the marker from its known pixel location(s) in the image. Alternatively, the transform-domain version of the known marker could be removed before the inverse transform of block  930 . Block  950  calculates the real part of the approximate restored image. Block  960  calculates the exact form of the restored image similarly to block  440 ,  FIG. 4 . Here again, the entire transformed version, including the magnitude component, of the restored image could be calculated if desired. 
       FIG. 10  shows a Matlab® program  1000  for restoring an image from a phase component, according to method  900 . After inverse-transforming this component, the program removes the marker phase phi and magnitude (radius) R from the base-domain form JC. After calculating the real part of the image J(1 . . . N 1 ,1 . . . N 2 ) from which the marker has been removed, the program applies the procedure described in connection with method  200  for producing the exact pixel values of the restored image in the base domain. 
     Restoration from Sign-Magnitude Component 
       FIG. 11  is a flow chart of a method  1100  for modifying a transform-domain image for communication as a single component in the transform domain. This example describes a one-dimensional row or column of an image X(1 . . . N)=ReX(1 . . . N)+jImX(1 . . . N), which may have multiple dimensions. The base-domain pixels are numbers that are rational in terms of a preselected precision, which may differ for different embodiments. 
     Block  1110  receives the image from an image-generating device, a storage, or other input. Normally, this is already expressed in the transform domain, as
 
   X   (1  . . . K )= Re X   (1  . . . K )+ jIm X   (1  . . . K )=   MX   (1  . . . K )exp( j PX   (1  . . . K )),
 
in Cartesian or polar coordinates. However, it could be received in the base domain and transformed in block  1110 ; the remainder of method  110 , and all of method  1300  below, would remain unaffected.
 
     Block  1120  generates a known marker having a strictly real value of magnitude  C . Appropriate values of  C  may be selected as described in connection with method  200 ,  FIG. 2 . 
     Block  1130  embeds marker  C  within image  X  to produce a modified original image
 
   XC   (1  . . . K )={ Re[ X   (1  . . . K )]+ C}+j{Im[ X   (1  . . . K )]}.
 
Because  C  is real, the imaginary part of  X  remains unchanged. The marker in this example represents the transform-domain equivalent of a delta functional in the base domain. As in method  200 , a base-domain marker C=cδ(n 1 ,n 2 ) could alternatively be applied to a pixel at location X(n 1 ,n 2 ) of the base-domain image, if that image is available at the source. For an image formed as a sequence of discrete pixels, applying a delta functional amounts to adding a predefined constant value to only one of the existing pixels, conveniently to the first element, such as X(1) in a one-dimensional image, or a row or column of a two-dimensional image. Again, markers may be functions other than transform-domain constants or base-domain delta functionals.
 
     Block  1140  calculates both the magnitude and phase of the modified transform-domain image from its real and imaginary parts. The magnitude component is
 
   MXC   (1  . . . K )=[{ Re[ XC   (1  . . . K )]})^2 +{Im[ XC   (1  . . . K )]}^2]^0.5.
 
The phase component is
 
   PXC   =arctan{ Im[ XC   (1  . . . K )]/ Re[ XC   (1  . . . K )]}.
 
     Block  1150  generates a sign map of  XC . This quantity is defined as  map (1 . . . K) =1−2×[ PXC (1 . . . K)&lt;0]. That is, each element of  map  has a value of +1 or −1 depending upon whether the corresponding element  PXC  has a phase that is negative or non-negative. Block  1160  embeds elements of the sign map into corresponding elements of the magnitude component,  SMXC (1 . . . K)= map * MXC , where the asterisk denotes term-by-term multiplication. 
     Block  1170  then sends only the magnitude component  SMXC  to medium  130 ,  FIG. 1 , for transmission and/or storage, and sends with it the sign map  map . Again, if the marker  C  is redetermined for each image, its value may be sent along with the signed-magnitude component. 
       FIG. 12  shows a Matlab® program for an example of a program  1200  for modifying an original image for communication as a single signed-magnitude component in a transform domain. The variable names are noted in the figure. 
       FIG. 13  illustrates a method  1300  for restoring an original base-domain image from only a signed-magnitude component. Block  1310  receives the signed magnitude component  SMXC (1 . . . K) from medium  130 ,  FIG. 1 . Block  1320  extracts the sign map from the sign of the component,  map (1 . . . .K)=1−2×[ SMXC (1 . . . K)&lt;0]. Block  1330  calculates the component&#39;s magnitude
     MXC =| SMXC |=[Re SMXC   (1  . . . K )^2 +Im SMXC   (1  . . . K )^2]^0.5. 
     Block  1340  expands the magnitude by one place of accuracy, in this example one decimal place, in the base domain of the image. Inverse transforming, EQA=ifft( MXC ). Then the base-domain part acquires an extra place, Integer( 2 ×EQA)/2=EQA. This quantity is then converted back to the transform domain, Re( XC )=fft(EQA). 
     Block  1350  restores the phase of the image. The raw phase is calculated as
 
   PXC   =arccos[ Re (   XC   )/   MXC   ].
 
Then the map supplies the actual phase,
 
   PSXC   (1  . . . K )=   map   (1  . . . K )*   PXC   (1  . . . K ).
 
Again, the asterisk denotes term-by-term multiplication.
 
     Block  1360  produces a base-domain image from the entire transform form, QCA=ifft[ MXC exp(j PXC )]. 
     Block  1370  removes the marker as in method  200 , QA(1)=QCA− C . When the marker is a delta functional, its value is subtracted only from the base-domain element where it was applied—the first element, in this example. 
     Block  1380  produces the exact base-domain image by rounding the approximate form, Q(1 . . . K)=Integer[QA(1 . . . K)]. This image has the exact element values of the original image in its base-domain form. For example, if the transform image received in block  1110  represented a 256-level gray-scale image in the base domain, then restored image Q has exactly the same accuracy. Block  1390  stores or transmits this restored image Q for output or further processing. 
       FIG. 14  shows a program  1400  for restoring the original image from a signed magnitude component. In this example, the original input image has two dimensions, and is received in the base domain, so that the original and restored images may be compared. 
     Complex Image Pixels 
     Previously described methods assumed that the pixels in the original image have only real values. Some applications, however, produce pixels that have complex values. Synthetic-aperture radars, seismic analysis, and some medical imaging, for example, produce transform-domain images that translate back to a spatial or other domain as complex-valued pixels. In addition, there are applications in which the ability to process complex pixel values is advantageous even when the actual pixel values are strictly real. 
       FIG. 15  is a flow chart of an example method  1500  for modifying such an image for communication as a single component in the transform domain, where the original image-pixel values may assume complex values. Block  1510  represents an original input image having complex pixel values in a base domain. Although the image may have any number of dimensions, this example focuses upon a two-dimensional matrix of pixels I(n 1 ,n 2 )=I 1 (n 1 ,n 2 )+j I 2 (n 1 ,n 2 ), representing the real and imaginary parts. Block  1511  transforms the base-domain image to a transform-domain image  X (k 1 ,k 2 )=fft 2 (I(n 1 ,n 2 ). If the original image arrives in transform-domain form, block  1511  is not needed. Unless otherwise from the context, equations are expressed in the syntax of the aforementioned Matlab® language. 
     Block  1520  generates a sign map similarly to block  720 ,  FIG. 7 ; because this example method employs other maps as well, this map is called a phase map. Block  1520  calculates the phase component  PX (k 1 ,k 2 ) of the transformed image and constructs a binary-valued map SGN 1 (k 1 ,k 2 ) having one value (e.g., “1”) where the phase has a negative value, and another value (e.g., “0”) where the sign of the phase is positive. 
     Block  1530  generates a real-valued marker, similarly to block  730 . An amplification factor may be employed, as mentioned previously, to ensure minimum phase of the modified image in the transform domain. This factor may be based upon the image pixel values, or may be a constant. Typical markers may be  C =α N 1 N 2 , where α=10^4 to 10^8. Block  1540  embeds the marker as described previously, and block  1550  calculates the magnitude component of the thus modified image,  MXC (k 1 ,k 2 )=| X (k 1 ,k 2 ) + C )|. 
     Blocks  1560  may be added to reduce image storage requirements in another way, by reducing the storage required for individual matrix values in the magnitude component. Block  1561  reduces the highest values of the matrix elements by subtracting the large marker value from every element,  MXC1 k 1 ,k 2 )= MXC (k 1 ,k 2 )− C . Because the resulting element values may be negative as well as positive, block  1561  then generates a sign map SGN 2 (k 1 ,k 2 ) having one value (e.g., “1”) for elements having a negative value—i.e., below a threshold value of 0—and another value (e.g., “0”) for positive elements, above 0 or some other threshold. Block  1562  reduces storage by allowing each matrix element of the magnitude component to occupy a variable number of bytes of storage. Block  1562  first computes the absolute values of the individual elements,  MXC1 (k 1 ,k 2 )=| MXC1 (k 1 ,k 2 ) |, and rounds them to the nearest integer. Block  1562  then generates a location map or byte-length map
 
BYTE( k   1   ,k   2 )=ceiling (|log 2(   MXC 1 ( k   1   ,k   2 ))/8|).
 
Each matrix element thus indicates how much storage is required for the corresponding element of  MXC1 . In this example, each element of BYTE has three bits, indicating that one through eight bytes may be used for each element of  MXC1 . This allows each of the elements of the magnitude matrix  MXC1 (k 1 ,k 2 ) to be stored with variable byte lengths, rather than requiring the same storage for each element, regardless of its actual magnitude. Storage units other than bytes may be employed here if desired. Although blocks  1560  add overhead, the storage/bandwidth savings may outweigh it in many cases.
 
     Block  1570  packages the three maps SGN 1 (k 1 ,k 2 ), SGN 2 (k 1 ,k 2 ), BYTE(k 1 ,k 2 ) along with the magnitude component  MXC1 (k 1 ,k 2 ). The maps may be embedded in the image itself, sent as separate files, or packaged in any other convenient manner. As before, the value of the marker  C  (or α. N 1 , N 2 ), as well as other information such as an image name, may be included if desired. Block  1580  represents the total package that is communicated—i.e., transmitted or stored for future restoration. 
       FIG. 16  depicts a method  1600  for restoring an original complex-valued base-domain image from a signed-magnitude transform-domain component  1610  such as that communicated in a medium  130 ,  FIG. 1 , from block  1580 ,  FIG. 15 . Block  1620  extracts the maps SGN 1 (k 1 ,k 2 ), SGN 2 (k 1 ,k 2 ), BYTE(k 1 ,k 2 ), and possibly other data, packaged along with th magnitude component  MXC 1. The phase (sign) and magnitude maps may be converted to a form that is more useful in this particular example, MAP 1 (k 1 ,k 2 )=1−2*SGN 1 (k 1 ,k 2 ) and MAP 2 (k 1 ,k 2 )=1−2*SGN 2 (k 1 ,k 2 ). Block  1630  reconstitutes the full magnitude component from the compressed form in which it was communicated to block  1610 . Location map BYTE(k 1 ,k 2 ) may optionally expand the variable-byte form by known methods to a form having the same number of bytes per element, for easier computation in later stages. The marker value subtracted in block  1561  is restored as MXC(k 1 ,k 2 )=MXC1(k 1 ,k 2 )+C. 
     Block  1640  calculates the real part of the restored transform-domain image. The marker is partially removed in the base domain during this process, as indicated at  1641 . The even part of the base-domain image is Jev(n 1 ,n 2 )=ifft 2 ( MXC ). An operation Jev(1,1)− C   removes the base-domain form of the real-valued marker, indicated as block  1641 . The real part of the transform-domain image is  RX =real (fft 2 (Jev)). 
     Block  1650  calculates the imaginary part of the restored transform-domain image—that is, restores the phase of the original image. Consider the following derivation of the imaginary part  IX  from  MXC , the modified magnitude: 
                     ⁢       MXC   ^   2     =             (       RX   _     +     C   _       )     ^   2     +       IX   _     ^   2       ⁢     
     ⁢           ∴           ⁢       IX   _     ^   2       =           MXC   _     ^   2     -       (       C   _     +     RX   _       )     ^   2       ⁢     
     ⁢           =         (       MXC   _     -     (       C   _     +     RX   _       )       )     *     (       MXC   _     +     (       C   _     +     RX   _       )       )       ⁢     
     ⁢           =         (       MXC   _     -     C   _     -     RX   _       )     *     (       MXC   _     +     C   _     +     RX   _       )       ⁢     
     ⁢           =           C   _     ^   2     *     (         MXC   _     /     C   _       -   1   +         R   ⁢           ⁢   X     _     /     C   _         )     *     
     ⁢           ⁢       (         MXC   _     /     C   _       +   1   +       RX   _     /     C   _         )     .     
     ⁢           ⁢   Defining     ⁢           ⁢     AN   _       ≡         MXC   _     /     C   _       -   1   +         RX   _     /     C   _       ⁢           ⁢   and                                       ⁢       AP   ≡         MXC   _     /     C   _       +   1   +       RX   _     /     C   _           ,     
     ⁢           ⁢       IX   _     =       C   _     *   sqrt   ⁢           ⁢     (     AN   _     )     *   sqrt   ⁢           ⁢     (     AP   _     )                 
This computation differs somewhat from earlier phase restoration in order to retain accuracy of the element values in view of precision limitations in the data processor, because  MXC  and  C  are very large numbers.  AN  should have all positive values. Block  1650  sets any slightly negative values due to precision shortcomings are to zero by  AN =( AN &gt;0).* AN . The correct (bipolar) version of  IX  is not needed in this example. If desired, the algebraic signs of the  IX  elements may be corrected by reintroducing the sign map, noted at block  1661 . The correction is  IX = IX .*RMAP.*MAP 1 , where RMAP=1−2*( RX &lt;0). The purpose of RMAP is to correct the signs that might have been changed due to precision errors.
 
     Block  1660  combines the real and imaginary parts to find the complete transform-domain restored image. In polar coordinates, the magnitude of the final transform-domain image then becomes  MX =sqrt( RX ^2+ IX ^2), and the phase is  PX =MAP 1 .* arccos( RX / MX ). Again, MAP 1  reintroduces the sign-map correction, indicated at block  1661 . If Cartesian coordinates  Y  are desired, the transform-domain image expression is  Y (k 1 ,k 2 )= RX (k 1 ,k 2 )+j  IX (k 1 ,k 2 ). 
     Block  1680  takes an inverse discrete Fourier transform to produce the base-domain restored image J(k 1 ,k 2 )=ifft 2 ( Y ). If only the transform-domain restored image is required, block  1680  may be omitted. 
     In methods  1500  and  1600 , the marker value  C  had a strictly real value. Modifications to those methods, however, would allow imaginary marker values, generated in the same way as the real values. In block  1550 , however, magnitude  MXC (k 1 ,k 2 ) would become  MXC (k 1 ,k 2 )=| X (k 1 ,k 2 )+j  C |. 
     In method  1600 , for an imaginary marker, the inverse transform of the magnitude matrix  MXC (k 1 ,k 2 ) produces an approximation to the odd portion Jod(n 1 ,n 2 ) of the complex image I(n 1 ,n 2 ), and the imaginary marker value  C  is subtracted from the appropriate element of the complex matrix. Where the marker resides in the first pixel, Jod(1,1)=Jod(1,1)−j  C . The imaginary part  IX (k 1 ,k 2 ) of the complex transform becomes  IX =−imag (fft 2 (Jod)). The matrix j  IX (k 1 ,k 2 ) is imaginary and bipolar, so a sign map may be calculated as IMAP(k 1 ,k 2 )=1−2 *( IX &lt;0). The real part  RX (k 1 ,k 2 ) of the of the complex transform matrix  X (k 1 ,k 2 ) remains the same. Retaining element precision here favors a slightly different computation:  BN = MXC / C −1+ IX / C . Again, any slightly negative values are set to zero by  BN =( BN &gt;0).* BN .  BN  corresponds to  AN , above.  BP = MXC / C +1+ IX / C . corresponds to  AP . The magnitude of the RX(k 1 ,k 2 ) may then be computed as RX=C*sqrt(BN).*sqrt(BP). The bipolar matrix version of the RX(k 1 ,k 2 ), if needed, is RX=IMAP.*MAP 1 .*RX. Block  1660  would calculate the phase component as PX(k 1 ,k 2 )=MAP 1 .*arcsin(IX/MX). Otherwise, the operations may be the same as described in connection with  FIGS. 15 and 16 . 
     Complex image elements, whether transform or base domain, are frequently expressed in single- or even double-precision floating-point notation from the beginning. Therefore in many cases an exact restoration down to the bit level is questionable. However, if the original image elements  X (k 1 ,k 2 ) are—or can be considered to be—integer values, then the rounding techniques described in connection with methods  400  and  800  may be employed in method  1600  as well. Although the storage-reduction blocks in methods  1500  and  1600  are especially useful in high-precision image processing, they may be incorporated in nay of the other methods as well. 
       FIG. 17  shows an application of methods capable of modifying and restoring complex-valued images. Method  1700  modifies, communicates, and restores a pair of independent real-valued images as a single component of a single image. Original images  1711  and  1712  may comprises images of any type, as long as their pixels or other elements do not have complex values. These will be termed “real-valued” pixels, although their elements could be strictly imaginary or represent any other single value. They are symbolized herein as I1(n 1 ,n 2 ) and I2(n 1 ,n 2 ). 
     Block  1720  converts one of the images to have imaginary pixel values, j*I2(n 1 ,n 2 ). Block  1730  then combines the images into a single image having complex pixel values, I(n 1 ,n 2 )=I1(n 1 ,n 2 )+j*I2(n 1 ,n 2 ). That is, each pixel of image  1711  become the real part of a pixel of the combined image, and each pixel of image  1712  becomes the complex part of a corresponding pixel in the combined image. These terms may include Cartesian angles other than 0° and 90°, or magnitude and phase, or any other pair of mutually orthogonal complex representations. 
     Block  1740  modifies the combined image I(n 1 ,n 2 ) to produce a single component (magnitude or phase)  1750 , which has the storage and bandwidth requirements of other modified images described above. This component may be communicated by transmission or storage on a medium such as 130,  FIG. 1 . 
     Block  1760  retrieves the single component and restores it according to methods such as 800 described above to produce the original combined image I(n 1 ,n 2 ). Block  1770  separates the restored image into its real and imaginary parts, ReI(n 1 ,n 2 ) and ImI(n 1 ,n 2 ). Block  1780  reconverts the imaginary part to real values, if necessary, say as I2(n 1 ,n 2 )=ImI(n 1 ,n 2 )/j. Blocks  1791  and  1792  represent the restored images I1(n 1 ,n 2 )=ReI(n 1 ,n 2 ), and I2(n 1 ,n 2 )=ImI(n 1 ,n 2 )/j. 
     CONCLUSION 
     Embodiments of the invention offer methods and systems for communicating images efficiently by reducing storage space and transmission time/bandwidth. Images may have any number of dimensions or formats, the term “image” herein is not limited to spatial images, and may include any sequence of data that is representable in an appropriate form. Complex or other multi-part element values may be accommodated. Any of the embodiments may process an original image that is presented either in a base domain or in a transform domain. Transforms other than DFT may be employed. The invention may be used along with other compression methods, and is transparent to them. 
     The foregoing description and drawing illustrate specific embodiments of the invention sufficiently to enable those skilled in the art to practice it. Other embodiments may incorporate structural, logical, electrical, process, and other changes. Examples merely typify possible variations. Individual components and functions are optional unless explicitly required, and the sequence of operations may vary, or may be performed in parallel. Portions and features of some embodiments may be included in or substituted for those of others. Items in lists are not exclusive, and may also be combined in any way. The Abstract is provided solely as a search aid, and is not to be used for claim interpretation. The scope of the invention encompasses the full ambit of the following claims and all available equivalents.