AU2008211991A1 - Frequency estimation under affine distortion - Google Patents

Description

S&F Ref: 862327 AUSTRALIA PATENTS ACT 1990 COMPLETE SPECIFICATION FOR A STANDARD PATENT Name and Address Canon Kabushiki Kaisha, of 30-2, Shimomaruko 3 of Applicant: chome, Ohta-ku, Tokyo, 146, Japan Actual Inventor(s): Kieran Gerard Larkin, Peter Alleine Fletcher Address for Service: Spruson & Ferguson St Martins Tower Level 35 31 Market Street Sydney NSW 2000 (CCN 3710000177) Invention Title: Frequency estimation under affine distortion The following statement is a full description of this invention, including the best method of performing it known to me/us: 5845c(1392311 _1) FREQUENCY ESTIMATION UNDER AFFINE DISTORTION FIELD OF THE INVENTION The present invention relates to spectral estimation and, in particular, the measurement of the Optical Transfer Function at designated frequencies that have undergone unknown affine distortions. 5 BACKGROUND Spectral estimation is a technique with very wide application. For example, in chemical spectroscopy the dominant spectral lines relate to certain atomic electron transitions, and the combination of lines represents a spectral "fingerprint" of the chemical to under test. In the analysis of sunspots and other long (or short) time series data with hidden periodicities, Fourier transform based spectral analysis is often used. In the case of equally spaced or sampled data, the dramatic speed of the fast Fourier transform (FFT) can often be used to great benefit in spectral estimation. One major drawback of the FFT based spectral estimators is known as the is windowing effect (or equivalently, spectral leakage). This effect is related to the fact that actual data is only available for finite length sequences, not necessarily an integer multiple of the underlying period. Sharp truncation of the dataset before applying Fourier methods introduces artefacts related to the sharp discontinuities. One approach to this effect is to apply a further, multiplicative, window to the dataset, but ensuring that the window has 20 smooth, gradual transition from fully on (1.0) to fully off (0.0). More sophisticated methods exist for ameliorating the spectral leakage effects of finite data lengths. In astronomy, and in particular radio astronomy, celestial images are derived from the two-dimensional (2-D) Fourier transform of signals from arrays of aerials. In effect, the desired image is the (2-D spatial) spectrum of a detected radio wave. Because 25 of restrictions on the size of the receiving aerial array, the spectrum-derived image is often severely distorted by spectral leakage effects, as well as other common effects like thermal and photon noise. One conventional approach is to apply strong constraints on the sparsity of the spectral image. For example, stars are expected to be essentially point-like. The method then works by subtracting the largest peak, scaled by the point spread function 13818451 DOC IRN: 862327 -2 (PSF) of the system, then moving on to the next peak and repeating the process until all peaks within some defined range have been removed, or until the residual meets some minimum criterion. There are a number of other methods (e.g., maximum entropy, Richardson-Lucy, and Wiener filtering to mention just three) with various advantages and 5 disadvantages. Two criteria by which a method can be judged are speed and accuracy. Algorithm complexity is sometimes important too, but is inversely related to the speed in most systems. Accuracy is important for photometric applications; where one needs to know the actual spectral power related to each peak. In general, it is not easy to maximise both speed and accuracy, because accuracy often requires many levels of iteration. 10 Thus, a need exists to provide an improved spectral estimation method. In particular, a need exists to provide a high-precision spectral estimation method. A further need exists to provide a high-speed spectral estimation method. SUMMARY 15 According to an aspect of the present disclosure, there is provided a measurement system and method for measuring the spectral amplitudes in systems with input images containing known spectral component geometries (or constellations), apart from unknown affine distortion. The proposed estimation process and system are able to determine at least one of the affine distortion and the spectral amplitude. 20 According to a first aspect of the present disclosure, there is provided a method of spectral estimation utilising a recorded image of a known test pattern, wherein the test pattern contains a known constellation of spectral peaks. The method includes the steps of: determining probable peaks in the recorded image; adding the probable peaks to a list of provisional peaks; performing an affine invariant search of all permutations of provisional 25 peaks in the list of provisional peaks using a least squares measure of best fit to identify a best fit; deriving affine parameters from the identified best fit; applying the affine parameters to an ideal constellation to generate a best fit constellation; determining using matrix algebra the interaction of all the peaks using a Fourier transform of a window function; disentangling peak values from the determined interactions; and utilising the 30 disentangled peak values to represent magnitudes and phases of a spectral estimate. 1381845 LDOC IRN 862327 -3 According to a second aspect of the present disclosure, there is provided a spectral estimation system for determining a spectral estimation of an imaging system by utilising a recorded image of a known test pattern, wherein said recorded image was recorded by said imaging system. The spectral estimation system includes an input for receiving said 5 recorded image, a storage medium for storing a provisional list of provisional peaks; and a processor. The processor performs the operations of: determining probable peaks in said recorded image; adding said probable peaks to a list of provisional peaks; performing an affine invariant search of all permutations of provisional peaks in said list of provisional peaks using a least squares measure of best fit to identify a best fit; deriving affine io parameters from said identified best fit; applying said affine parameters to an ideal constellation to generate a best fit constellation; determining using matrix algebra interaction of all the peaks using a Fourier transform of a window function; and disentangling peak values from said determined interactions. The spectral estimation system further includes an output for presenting said disentangled peak values as 15 representations of magnitudes and phases of a spectral estimate. According to a third aspect of the present disclosure, there is provided a computer readable storage medium having a computer program recorded thereon, the program being executable by a computer apparatus to make the computer perform a method of spectral estimation utilising a recorded image of a known test pattern, wherein said test pattern 20 contains a known constellation of spectral peaks. The program includes: code for determining probable peaks in said recorded image; code for adding said probable peaks to a list of provisional peaks; code for performing an affine invariant search of all permutations of provisional peaks in said list of provisional peaks using a least squares measure of best fit to identify a best fit; code for deriving affine parameters from said 25 identified best fit; code for applying said affine parameters to an ideal constellation to generate a best fit constellation; code for determining using matrix algebra interaction of all the peaks using a Fourier transform of a window function; code for disentangling peak values from said determined interactions; and code for utilising the disentangled peak values to represent magnitudes and phases of a spectral estimate. 30 1381845 .DOC IRN: 862327 -4 BRIEF DESCRIPTION OF THE DRAWINGS One or more embodiments of the invention will now be described with reference to the drawings, in which: Fig. 1 shows a plot of an intensity function of a typical spatial greyscale test pattern; 5 Fig. 2 shows a plot of the Fourier spectrum of the test pattern of Fig. 1; Fig. 3 shows a plot of a Fourier spectrum of an ideal test pattern; Fig. 4 shows a plot of a Fourier spectrum corresponding to a square windowed test pattern; Figs 5A - C show plots of the real, imaginary and magnitude spectral components 1o for a test pattern with two (interacting) sinusoids; Figs 6A - C show a selection of some possible window functions; Fig. 7 is a flow diagram of a complete spectral estimation process in accordance with an embodiment of the present disclosure; and Figs 8A and 8B form a schematic block diagram of a general purpose computer is system upon which the arrangements described can be practised. DETAILED DESCRIPTION Where reference is made in any one or more of the accompanying drawings to steps and/or features, which have the same reference numerals, those steps and/or features have 20 for the purposes of this description the same function(s) or operation(s), unless the contrary intention appears. Overview One method in accordance with an embodiment of the present disclosure is restricted to a special class of spectral models with pre-ordained geometry. Constraining the spectral 25 model is not possible in more general applications, such as astronomy. In such scenarios, it is unreasonable to expect that the spatial arrangement of stars be defined before measurement. Nevertheless, in many spectral estimation situations it may be possible to pre-ordain certain aspects of the spectrum. One such circumstance is the spectral estimation of imaging systems in which predefined test charts are used. Typical examples 30 of such imaging systems are digital cameras, document scanners, and the lenses of digital 1381845 _1 I IRN: 862327 -5 SLR cameras. Spectral estimation allows the measurement of the input-output contrast at various spatial frequencies as a way of characterising the spatial frequency response (resolution) of the system. The nominal geometry can undergo an unknown affine distortion without adversely affecting the estimation process, if the estimation process is 5 effectively affine invariant. According to an aspect of the present disclosure, there is provided a measurement system and method for rapidly and precisely measuring the spectral amplitudes in systems with input images containing known spectral component geometries (or constellations), apart from unknown affine distortion. The proposed estimation process and system find 10 both the affine distortion and the spectral amplitude. According to one embodiment of the present disclosure, there is provided a method of spectral estimation utilising a recorded image of a known test pattern, wherein the test pattern contains a known constellation of spectral peaks. The method includes the steps of: determining probable peaks in the recorded image; adding the probable peaks to a list of is provisional peaks; performing an affine invariant search of all permutations of provisional peaks in the list of provisional peaks using a least squares measure of best fit to identify a best fit; deriving affine parameters from the identified best fit; applying the affine parameters to an ideal constellation to generate a best fit constellation; determining using matrix algebra the interaction of all the peaks using the Fourier transform of the window 20 function; disentangling peak values from the determined interactions; and utilising the disentangled peak values to represent magnitudes and phases of a spectral estimate. According to another embodiment of the present disclosure, there is provided a spectral estimation system for determining a spectral estimation of an imaging system by utilising a recorded image of a known test pattern, wherein said recorded image was 25 recorded by said imaging system. The spectral estimation system includes an input for receiving said recorded image, a storage medium for storing a provisional list of provisional peaks; and a processor. The processor performs the operations of: determining probable peaks in said recorded image; adding said probable peaks to a list of provisional peaks; performing an affine invariant search of all permutations of provisional peaks in 30 said list of provisional peaks using a least squares measure of best fit to identify a best fit; deriving affine parameters from said identified best fit; applying said affine parameters to 1381845 1.DOC IRN: 862327 -6 an ideal constellation to generate a best fit constellation; determining using matrix algebra the interaction of all the peaks using the Fourier transform of the window function; and disentangling peak values from said determined interactions. The spectral estimation system further includes an output for presenting said disentangled peak values as 5 representations of magnitudes and phases of a spectral estimate. According to another embodiment of the present disclosure, there is provided a computer readable storage medium having a computer program recorded thereon, the program being executable by a computer apparatus to make the computer perform a method of spectral estimation utilising a recorded image of a known test pattern, wherein 10 said test pattern contains a known constellation of spectral peaks. The program includes: code for determining probable peaks in said recorded image; code for adding said probable peaks to a list of provisional peaks; code for performing an affine invariant search of all permutations of provisional peaks in said list of provisional peaks using a least squares measure of best fit to identify a best fit; code for deriving affine parameters from said is identified best fit; code for applying said affine parameters to an ideal constellation to generate a best fit constellation; code for determining using matrix algebra interaction of all the peaks using a Fourier transform of a window function; code for disentangling peak values from said determined interactions; and code for utilising the disentangled peak values to represent magnitudes and phases of a spectral estimate. 20 In one embodiment, the constellation defines the number and location of the spectral peaks. In a further embodiment, the recorded image was recorded by an imaging system. In one implementation, the imaging system is chosen from the group of imaging systems including digital cameras, document scanners, and the lenses of digital SLR cameras. In one emobiment, the spectal estimate characterises at least one property of the imaging 25 system used to record the test pattern. 1381845 .DOC IRN: 862327 -7 Computer Implementation Figs. 8A and 8B collectively form a schematic block diagram of a general purpose computer system 800, also referred to as a computer apparatus, upon which the various arrangements described can be practised. 5 As seen in Fig. 8A, the computer system 800 is formed by a computer module 801, input devices such as a keyboard 802, a mouse pointer device 803, a scanner 826, a camera 827, and a microphone 880, and output devices including a printer 815, a display device 814 and loudspeakers 817. An external Modulator-Demodulator (Modem) transceiver device 816 may be used by the computer module 801 for communicating to and 10 from a communications network 820 via a connection 821. The network 820 may be a wide-area network (WAN), such as the Internet or a private WAN. Where the connection 821 is a telephone line, the modem 816 may be a traditional "dial-up" modem. Alternatively, where the connection 821 is a high capacity (eg: cable) connection, the modem 816 may be a broadband modem. A wireless modem may also be used for wireless 15 connection to the network 820. The computer module 801 typically includes at least one processor unit 805, and a memory unit 806 for example formed from semiconductor random access memory (RAM) and semiconductor read only memory (ROM). The module 801 also includes an number of input/output (1/0) interfaces including an audio-video interface 807 that couples to the 20 video display 814, loudspeakers 817 and microphone 880, an I/O interface 813 for the keyboard 802, mouse 803, scanner 826, camera 827 and optionally a joystick (not illustrated), and an interface 808 for the external modem 816 and printer 815. In some implementations, the modem 816 maybe incorporated within the computer module 801, for example within the interface 808. The computer module 801 also has a local network 25 interface 811 which, via a connection 823, permits coupling of the computer system 800 to a local computer network 822, known as a Local Area Network (LAN). As also illustrated, the local network 822 may also couple to the wide network 820 via a connection 824, which would typically include a so-called "firewall" device or device of similar functionality. The interface 811 may be formed by an EthernetTM circuit card, a 30 BluetoothTM wireless arrangement or an IEEE 802.11 wireless arrangement. 1381845 LDOC IRN: 862327 -8 The interfaces 808 and 813 may afford either or both of serial and parallel connectivity, the former typically being implemented according to the Universal Serial Bus (USB) standards and having corresponding USB connectors (not illustrated). Storage devices 809 are provided and typically include a hard disk drive (HDD) 810. Other storage 5 devices such as a floppy disk drive and a magnetic tape drive (not illustrated) may also be used. An optical disk drive 812 is typically provided to act as a non-volatile source of data. Portable memory devices, such as optical disks (eg: CD-ROM, DVD), USB-RAM, and floppy disks, for example, may then be used as appropriate sources of data to the system 800. to The components 805 to 813 of the computer module 801 typically communicate via an interconnected bus 804 and in a manner which results in a conventional mode of operation of the computer system 800 known to those in the relevant art. Examples of computers on which the described arrangements can be practised include IBM-PCs and compatibles, Sun Sparcstations, Apple MacTM or alike computer systems evolved 15 therefrom. The method of spectral estimation may be implemented using the computer system 800 wherein the processes of Figs 1 to 7, to be described, may be implemented as one or more software application programs 833 executable within the computer system 800. In particular, the steps of the method of spectral estimation are effected by 20 instructions 831 in the software 833 that are carried out within the computer system 800. The software instructions 831 may be formed as one or more code modules, each for performing one or more particular tasks. The software may also be divided into two separate parts, in which a first part and the corresponding code modules perform the spectral estimation methods and a second part and the corresponding code modules manage 25 a user interface between the first part and the user. The software 833 is generally loaded into the computer system 800 from a computer readable medium, and is then typically stored in the HDD 810, as illustrated in Fig. 8A, or the memory 806, after which the software 833 can be executed by the computer system 800. In some instances, the application programs 833 may be supplied to the user 30 encoded on one or more CD-ROM 825 and read via the corresponding drive 812 prior to storage in the memory 810 or 806. Alternatively the software 833 may be read by the 1381845 .DOC RN: 862327 -9 computer system 800 from the networks 820 or 822 or loaded into the computer system 800 from other computer readable media. Computer readable storage media refers to any storage medium that participates in providing instructions and/or data to the computer system 800 for execution and/or s processing. Examples of such storage media include floppy disks, magnetic tape, CD-ROM, a hard disk drive, a ROM or integrated circuit, USB memory, a magneto-optical disk, or a computer readable card such as a PCMCIA card and the like, whether or not such devices are internal or external of the computer module 801. Examples of computer readable transmission media that may also participate in the provision of software, 10 application programs, instructions and/or data to the computer module 801 include radio or infra-red transmission channels as well as a network connection to another computer or networked device, and the Internet or Intranets including e-mail transmissions and information recorded on Websites and the like. The second part of the application programs 833 and the corresponding code modules is mentioned above may be executed to implement one or more graphical user interfaces (GUIs) to be rendered or otherwise represented upon the display 814. Through manipulation of typically the keyboard 802 and the mouse 803, a user of the computer system 800 and the application may manipulate the interface in a functionally adaptable manner to provide controlling commands and/or input to the applications associated with 20 the GUI(s). Other forms of functionally adaptable user interfaces may also be implemented, such as an audio interface utilizing speech prompts output via the loudspeakers 817 and user voice commands input via the microphone 880. Fig. 8B is a detailed schematic block diagram of the processor 805 and a "memory" 834. The memory 834 represents a logical aggregation of all the memory 25 devices (including the HDD 810 and semiconductor memory 806) that can be accessed by the computer module 801 in Fig. 8A. When the computer module 801 is initially powered up, a power-on self-test (POST) program 850 executes. The POST program 850 is typically stored in a ROM 849 of the semiconductor memory 806. A program permanently stored in a hardware device such as 30 the ROM 849 is sometimes referred to as firmware. The POST program 850 examines hardware within the computer module 801 to ensure proper functioning, and typically 1381845 LDOC IRN: 862327 -10 checks the processor 805, the memory (809, 806), and a basic input-output systems software (BIOS) module 851, also typically stored in the ROM 849, for correct operation. Once the POST program 850 has run successfully, the BIOS 851 activates the hard disk drive 810. Activation of the hard disk drive 810 causes a bootstrap loader program 852 5 that is resident on the hard disk drive 810 to execute via the processor 805. This loads an operating system 853 into the RAM memory 806 upon which the operating system 853 commences operation. The operating system 853 is a system level application, executable by the processor 805, to fulfil various high level functions, including processor management, memory management, device management, storage management, software 10 application interface, and generic user interface. The operating system 853 manages the memory (809, 806) in order to ensure that each process or application running on the computer module 801 has sufficient memory in which to execute without colliding with memory allocated to another process. Furthermore, the different types of memory available in the system 800 must be used is properly so that each process can run effectively. Accordingly, the aggregated memory 834 is not intended to illustrate how particular segments of memory are allocated (unless otherwise stated), but rather to provide a general view of the memory accessible by the computer system 800 and how such is used. The processor 805 includes a number of functional modules including a control 20 unit 839, an arithmetic logic unit (ALU) 840, and a local or internal memory 848, sometimes called a cache memory. The cache memory 848 typically includes a number of storage registers 844 - 846 in a register section. One or more internal buses 841 functionally interconnect these functional modules. The processor 805 typically also has one or more interfaces 842 for communicating with external devices via the system 25 bus 804, using a connection 818. The application program 833 includes a sequence of instructions 831 that may include conditional branch and loop instructions. The program 833 may also include data 832 which is used in execution of the program 833. The instructions 831 and the data 832 are stored in memory locations 828-830 and 835-837 respectively. Depending 30 upon the relative size of the instructions 831 and the memory locations 828-830, a particular instruction may be stored in a single memory location as depicted by the 1381845 .DOC IRN: 862327 - 11 instruction shown in the memory location 830. Alternately, an instruction may be segmented into a number of parts each of which is stored in a separate memory location, as depicted by the instruction segments shown in the memory locations 828-829. In general, the processor 805 is given a set of instructions which are executed therein. 5 The processor 805 then waits for a subsequent input, to which it reacts to by executing another set of instructions. Each input may be provided from one or more of a number of sources, including data generated by one or more of the input devices 802, 803, data received from an external source across one of the networks 820, 822, data retrieved from one of the storage devices 806, 809 or data retrieved from a storage medium 825 inserted 1o into the corresponding reader 812. The execution of a set of the instructions may in some cases result in output of data. Execution may also involve storing data or variables to the memory 834. The disclosed spectral estimation arrangements use input variables 854, that are stored in the memory 834 in corresponding memory locations 855-858. The spectral 15 estimation arrangements produce output variables 861, that are stored in the memory 834 in corresponding memory locations 862-865. Intermediate variables may be stored in memory locations 859, 860, 866 and 867. The register section 844-846, the arithmetic logic unit (ALU) 840, and the control unit 839 of the processor 805 work together to perform sequences of micro-operations 20 needed to perform "fetch, decode, and execute" cycles for every instruction in the instruction set making up the program 833. Each fetch, decode, and execute cycle comprises: (a) a fetch operation, which fetches or reads an instruction 831 from a memory location 828; 25 (b) a decode operation in which the control unit 839 determines which instruction has been fetched; and (c) an execute operation in which the control unit 839 and/or the ALU 840 execute the instruction. Thereafter, a further fetch, decode, and execute cycle for the next instruction may be 30 executed. Similarly, a store cycle may be performed by which the control unit 839 stores or writes a value to a memory location 832. 1381845_1 DOC IRN: 862327 - 12 Each step or sub-process in the processes of Figs I to 7 is associated with one or more segments of the program 833, and is performed by the register section 844-847, the ALU 840, and the control unit 839 in the processor 805 working together to perform the fetch, decode, and execute cycles for every instruction in the instruction set for the noted 5 segments of the program 833. The method of spectral estimation may alternatively be implemented in dedicated hardware such as one or more integrated circuits performing the functions or sub functions of determining probable peaks, adding probable peaks to a list of provisional peaks, performing an affine invariant search of all permutations of provisional peals to identify a 10 best fit, deriving affine parameters from the best fit, applying the affine parameters to an ideal constellation to generate a best fit constellation, determining the interaction of all peaks using a Fourier transform of the window function, disentangling peak values, and utilising the disentangled peak values to represent magnitudes and phases of a spectral estimate. Such dedicated hardware may include graphic processors, digital signal is processors, or one or more microprocessors and associated memories. Test Pattern The starting point for the spectral estimation problem is a pre-defined spectral test pattern. The test pattern is composed of a plurality of discrete frequency components. In 20 the case of a two-dimensional test pattern, at least two separate frequency components are necessary to define uniquely an affine transformation (excluding translation). The two components must not be collinear with the DC, where DC is an abbreviation for Direct Current or zero frequency component. In practice, the number of frequency components, denoted here by K, may be considerably greater than two. Typically, the image is defined 25 on an N by N pixel grid (containing N2 pixels in total). The discrete frequency components may contain exact integer periods over the 2D array, or can in general have non-integer (i.e., fractional) periods. Note that after affine distortion, even integer frequencies can become non-integer frequencies. 1381845_1 DOC IRN: 862327 - 13 Spectrum (Fourier transform of) Test Pattern In the ideal, continuous case, in which a test pattern extends infinitely, the Fourier transform of the test pattern is a collection of 2K delta functions. The doubling in the number is due to the fact that all real functions have Fourier transforms with Hermitian 5 symmetry and so (except for the special case of zero frequency or DC) any frequency component must have a Hermitian twin. Fig. I shows a plot 100 of an intensity function 150 of a 2-D test pattern. The plot 100 has an Intensity axis 110, a Spatial Dimension X axis 130 and a Spatial Dimension Y axis 140. In the plot 100, the Intensity axis 110 has a range of 0 to 6, the 10 Spatial Dimension X axis 130 has a range of 0 to I and the Spatial Dimension Y axis 140 has a range of 0 to 1. The plot 100 provides a texture map of the test pattern. Fig. 2 shows a plot 200 that shows the magnitude of the spectrum (namely the Fourier transform of the test pattern) of the intensity function 150 of Fig. 1. The plot 200 has an Spectral Magnitude axis 205, a Frequency Dimension u axis 215 and a is Frequency Dimension v axis 225. In the plot 200, the Spectral Magnitude axis 205 has a range of 0 to 0.8, the Frequency Dimension u axis 215 has a range of -10 to 10 and the Frequency Dimension v axis 225 has a range of -10 to 10. The plot 200 includes 4 symmetric pairs of peaks 210, 220, 230, and 240 corresponding to the 4 sinusoidal components of the texture shown in Fig. 1. 20 In an actual system, the test pattern is finite, and is only sampled on a discrete grid. The finite extent of the test pattern gives rise to extended functions (rather than delta functions) in the Fourier spectrum. If the test pattern is sharply limited (or windowed) to make it finite, the extended function will be of a since type. Other windows (as is well known in the literature) can have gentler cut-off or limiting, but result in a wider central 25 peak width in the Fourier transform. Alternatively, an even sharper window can have a narrower central peak, but the sidelobes then have greater power and oscillation. Windowing is a subject that has been extensively studied over the last half century or more. In general, windowing requires many compromises to be made. One embodiment of the present disclosure utilises all the information available equally over the entire window. 30 Furthermore, the embodiment utilises the extra information available, namely that there are a limited number of discrete spectral components. The problem then becomes one of 1381845 IDOC IRN: 862327 -14 precisely estimating the (possibly fractional) location of each spectral component and estimating the (typically complex) amplitude of the spectral component. In general, it is difficult to estimate the amplitude of a component without knowing its exact location. Conversely, it is difficult to estimate the location of a component 5 without knowing its exact amplitude. This is even true for a single component K=1, which gives rise to a pair of Hermitian conjugate components in the Fourier transform domain, 2K=2. Even if the exact location is known, the amplitude at any location is the sum of the overlapping sinc functions. In order to extract the single component amplitude, all the other interacting and overlapping amplitudes must be removed, and so their exact locations 1o must be known as well. In general, it is difficult to break this circular argument. An embodiment of the present disclosure utilises the prior (partial) knowledge of the spatial arrangement of spectral components as a starting point and builds from there. Fig. 3 shows a plot 300 of a Fourier sidelobe structure for an ideal test pattern. The 15 plot 300 has an Spectral Magnitude axis 305, a Frequency Dimension u axis 315 and a Frequency Dimension v axis 325. In the plot 300, the Spectral Magnitude axis 305 has a range of 0 to 0.75, the Frequency Dimension u axis 315 has a range of -6 to 0 and the Frequency Dimension v axis 325 has a range of 0 to 10. The plot 300 includes three clustered peaks 340, 350, 360, a separate peak 320 and a substantially flat surface 330. 20 Fig. 4 shows a plot 400 of a Fourier sidelobe structure for a square windowed test pattern. The plot 400 has an Spectral Magnitude axis 405, a Frequency Dimension u axis 415 and a Frequency Dimension v axis 425. In the plot 400, the Spectral Magnitude axis 405 has a range of 0 to 1, the Frequency Dimension u axis 415 has a range of-6 to 0 and the Frequency Dimension v axis 425 has a range of 0 to 10. The plot 400 includes 25 three clustered peaks 440, 450, 460, a separate peak 420 and a textured surface 430. Figs 5A, 5B, and 5C show plots of the real, imaginary, and magnitude components of two interacting sidelobes for a test pattern with two interacting sinusoids. Fig. 5A is a plot 510 of the real component, Fig. 5B is a plot 520 of the imaginary component, and Fig. 5C is a plot 530 of the magnitude component. 30 Fig. 7 is a flow diagram of a complete spectral estimation process 700 in accordance with an embodiment of the present disclosure. The spectral estimation process 700 begins 1381845 _1DOC IRN: 862327 - 15 at an initial step 710, which captures (records) a digital image of a known test pattern, using an imaging system that is to be characterised. As described above, typical examples of such imaging systems include digital cameras, document scanners, and the lenses of digital SLR cameras. In one embodiment, the test pattern contains a known constellation 5 of spectral peaks. Control passes to step 720, which selects a tile from the image captured in step 710. The tile is typically square and of size N byN pixels. However, the tile can equally be another shape, such as, for example, rectangular or hexagonal. Step 722 removes the DC component from the tile and control passes to step 724 to apply a chosen window function to the tile and zero pad the tile to a desired size. Two times padding in 10 each dimension is usually adequate. Control passes to step 726 to apply a 2-D discrete Fourier transform or fast Fourier transform (DFT or FFT). Step 728 detects L provisional peaks in the Fourier half-plane, thus determining probable peaks in the recorded image. It is not necessary to find Hermitian twins. Typically, L >> K. In one implementation, the provisional peaks are is added to a stored list of provisional peaks. The next steps of the process perform an affine invariant search of all permutations of provisional peaks in the list of provisional peaks using a least squares measure of best fit to identify a best fit and then derive affine parameters from the identified best fit. Thus, step 730 performs constellation matching by applying a least squares (LSF) fit to possible 20 affine transforms of ideal peak configuration. The constellation matching is repeated for all permutation of K samples chosen from L samples. Step 732 selects an affine permutation with the lowest residual error to be the best contender. Step 734 applies chosen affine parameters to ideal constellation positions of an ideal constellation, to generate a best fit constellation. 25 The method then uses matrix algebra to determine the interaction of all the peaks using the Fourier transform of the window function. The peak values are then disentangled from the determined interactions and the disentangled peak values are subsequently utilised 1381845 .DOC IRN: 862327 -16 to represent magnitudes and phases of a spectral estimate. Thus, step 736 generates a matrix for interaction of complex peaks located at positions defined by step 734. Control passes to step 738 to solve the matrix equation for interactions. This entails putting peak amplitudes (complex) into a linear matrix inversion equation defined by step 736. Step 5 740 uses the actual amplitudes and 2K component pairs to calculate a full 2D (NxN) amplitude pattern defined by discrete component. Step 742 subtracts the calculated pattern from a measured pattern to find the residual. Decision step 744 compares the residual with a reasonably expected residual to see if result is credible. If the residual is acceptable, Yes, a certain measure of confidence in the process is indicated and control proceeds to step to 746. If the residual is not acceptable, No, control proceeds to step 750. Decision step 746 determines whether the DC component of the residual is smaller than a preset threshold value. If Yes, control proceeds to a terminating Finish step 748. If the DC component of the residual is not smaller that the preset threshold value, No, control return to step 722 to remove the extra DC component before continuing on to step 724. is Decision step 750 determines whether to try a lower threshold. If Yes, control passes to step 728 to lower the threshold to detect whether more suitable peaks can be found. If so, use additional peaks in step 730 and continue process. If a lower threshold does not generate more useful peaks and is not to be tried at step 750, then control passes to a terminating Finish step 760. 20 At the end of the complete spectral estimation process 700 of Fig. 7, all the components, including DC, have been estimated and the residual image is known. The residual may contain noise as well as other artefacts, such as nonlinearities (harmonics and inter-modulation distortion). In one embodiment, the residual normalised by the total energies of the both the test 25 pattern image and the recorded image indicate a measure of the confidence in the process. The complete spectral estimation process 700 described above depends on the assumption that the spectral component configuration (constellation) is possible to estimate accurately regardless of the interactions. Note that interactions generally have two effects. 1381845 .DOC IRN: 862327 - 17 Each component peak location is shifted, and the resulting peak amplitude is also changed by the interaction. It is only the peak shift that causes problems. The amplitudes are automatically corrected by the linear matrix inversion step (steps 736 to 738), but only if the peak locations are correct. s There is an underlying assumption that a constellation consisting of a large number of spectral peaks can be located more accurately than a single spectral peak, or even a small number of spectral peaks. This will only be true if the peak shifts due to interaction tend to be (incoherently) in all directions rather than all coherently in one direction. It is also to be expected for constellations that contain components that are separated by many periods. A 10 fractional scaling or a rotation will tend to produce the desirable incoherent peak shifts. Pre-emptive DC Removal In most practical systems with real, positive images, the DC component is far larger than any other component, and exercises an excessive influence on all the other 15 components. To avoid some interactions in later process steps, the spectral estimation process of the present disclosure removes the DC level (or average) from the image tile as the first step, as described above with reference to step 722 of Fig. 7. The resulting image is then a zero-mean image. Later correction to the DC removal is possible, as described above with reference to step 746 of Fig. 7. 20 Fourier Analysis of Windowing The continuous function mathematical analysis will be described in this section for convenience and brevity, but note that all operations have their discrete data sample versions which map directly. The continuous Fourier Transform (FT) is in practice 25 approximated by the Discrete Fourier Transforms (DFTs) and will be usually evaluated using Fast Fourier Transforms (FFTs) because of the significant speed increases possible with the FFT. 1381845 .DOC IRN: 862327 -18 A first step defines an initial N-point pattern constellation by its Euclidean position vector components in Fourier space: (u,, ), n =1I -> N ... (1) Back in the spatial domain, this corresponds to an infinite texture test pattern 5 consisting of sine and cosine component: N f(x,y)= a, + (c. cos(21r[u,x+vyD+s, sin(21r[u,x+v,y} ... (2) The infinite pattern is actually only known or sampled over a finite domain (the tile) and the finite tile is further subjected to a window function w(x,y). Looking at the ideal test pattern in the Fourier domain, where the (continuous) Fourier transform is defined as 10 follows: F(u,v)= f ff(x,y)exp(2rd[ux+vyDdxdy ... (3) or in convenient mathematical shorthand: f (x,Y) 4 "T : F(u, V) ... (4) gives a discrete set of sampling Dirac delta functionals: F(u, v)= a 0 c(u, v)+ (,S(u -u,, v - v.)+5(u +u,,, v+ vJ) s2 (5) + {S(u - un, v - v.) - 5(u + u,, v + vA) 2i Transforming the windowed test pattern gives a rather different result. The well known Fourier convolution theorem is utilised, where * * represents 2-D convolution: g(x,y)= w(x, y).f(x, y)< " : G(u, v) = W(u, v)* *F(u, v) ... (6) 1381845_1 DOC IRN: 862327 -19 This leads on to the final form of the windowed tile FT: F(u,v)= aW(u,v)+ ( IW(u -u,,v- v.)+ W(u +u,,,v+ v.)) 27 2 ... (7) + L {W(u - u,, v - v") - W(u + u,, v + v.)} 2i which can be seen to be 2N + I copies of the window transform arranged on the constellation points. Because the window is finite and no bigger than the tile boundary, the s window transform has broad extent in Fourier space, and all the copies above interact additively with each other. Constellation Matching An aspect of the present disclosure is directed to determining and matching the io constellation points regardless of any affine geometric distortions. An advantage of constellation matching in the Fourier domain is that it is not necessary to consider shift (or translation) of the constellation, because a shift in the spatial domain does not cause a shift in the Fourier domain. In fact, a shift in the spatial domain only causes a phase change in the Fourier domain. This means that a six parameter affine distortion only requires a four 15 parameter affine search in the Fourier domain, and therefore two non-zero, non collinear, frequencies are the minimum necessary for affine determination. The relation between spatial and Fourier affine distortions is well known in the literature. The Fourier affine transform is defined by the 2x2 matrix: u' a b u, u, =au,+bv, (8) v, c d) v, v, = cu, + dv, 20 There are a number of methods for estimating the most likely correspondence between a set of ideal peaks in a certain constellation and set of detected peaks in an affine distortion of the ideal constellation. For example, it is possible to use some special properties of the ideal constellation, such as the arrangement of peaks along straight lines. Peak locations on straight lines are easily detected. However, disclosed herein is a more 1381845 I.DOC IRN: 862327 -20 general scheme that will work for arbitrary constellation geometries, and so can be utilised in all scenarios. Consider an ideal test pattern texture consisting of N real, non-zero sinusoidal components (and hence 2N +1 Fourier peaks in constellation). The method requires that s L > N provisional peaks are selected from the Fourier transform magnitude of the windowed tile. A subset N of the provisional peaks are selected and all permutations [total P(L, N) = LV(L - N)] of constellation matches with the N ideal peaks are evaluated. The permutation with the best match is chosen to define the affine distortion matrix. The best affine match criterion ((a,b,c,d) is evaluated as follows: to 4(a, b, c, d)= I(, -u')2 +(iv, -v')} ... (9) where the permuted constellation points are (i.,ii,) and the affine distorted ideals are (u',,v). The criterion has to be minimised in each instance (permutation) with respect to the four affine coefficients (a,b,c,d): 2 =0 ... (10) aa a N V " N = (, ,= ( -au,- ... (12) 8b ab 2 v(,, -v,) = (ii -cu,- dv.,, ,= 0 ... (12) ac ac 2 = (ii.-v',) =2 (ii.-cu,-dv.).v, =0 ... (13) ad ad 1381845 1.DOC IRN: 862327 - 21 These four equations determine the four affine unknowns and can be solved as follows: v2CUIi(= J,, - = iJuv (v14 a ( = 1 " "='2... (14) b = = " "- " "'(15) N NN (uj(v2 -(gunv, aboevle o t( q e ,,- uv) s0 cau= "=' ")"= " i 'c~~2 }-.. (18) n= = = )( u r ( , - ,( nv n d = "' "- " " .. ( 7 N NN E u " v - Eu,,v,, Now, the minimised total error for this affine permutation can be computed for the above values of the quartet (a, b, c,d): 10 E 2 -auln - bvn) + (V,, - CU, - dv} ... ( 18) Note that the above residual error can also be calculated directly from the variance covariance matrix. The process is repeated for all permutations and the permutation with the smallest match criterion is chosen to be the solution. Many variations on the above affine fitting process are also possible. For example, weighted least squares and Random 1381845_1.DOC IRN: 862327 - 22 Sample Consensus (RANSAC) can be utilised, and may be appropriate for certain special cases of templates and distortion processes, or larger numbers of peaks. At this stage of the process, the most likely affine distortion parameters and the most likely permutation of the constellation peaks are known. This means that the least squares 5 estimated peak locations (u,,, v,), and the corresponding (complex) peak values F(u,,v',) are known. These peak values are likely to be located at non-integer sample positions, and therefore should be estimated by high-order interpolation (such as Fourier interpolation) to ensure high accuracy. 1o Peak Detection The overall method is dependent on selecting a sufficient set of provisional peaks. Peaks are defined to be local maxima in the magnitude or more correctly the magnitude squared (or intensity) of the complex spectrum. Detecting M peaks is not the same as detecting the M largest values in the spectrum. A peak is defined by a zero gradient of the is intensity. A peak is surrounded on all sides by lower values of the intensity. Consequently, a practical peak finding method must define the size of region surrounding a prospective peak wherein all values are lower than the peak. This is essentially equivalent to defining a local scale parameter for peak occurrence. Peaks cannot then be located closer than this scale to adjacent peaks. Another way to view this is as a limit on the 20 maximum peak density. In practice, the task of peak detection can be converted to the more straightforward task of spiral phase detection. In one implementation, peak detection is implemented as spiral phase detection on the complex derivative of one of the spectral magnitude and intensity. Applying a complex derivative (Cauchy) operator to the spectral intensity 25 produces a complex result, the phase of which contains positive and negative spiral singularities. One polarity represents minima, the opposite polarity maxima. Spirals are then simply detected by circuit integrals around each point or pixel. The chose scale is simply determined by the size and shape of the circuit. 1381845 i.DOC IRN: 862327 -23 Compensating for the Window As shown in the windowing Fourier analysis, windowing causes the spectral peaks to spread in a complicated, but deterministic, manner. The spread function for each spectral peak is computed from the known window w(x,y). Typically, this is just a square array of 5 sample values. The most straightforward case involves all window sample values being 1.0, resulting in a Fourier transform W(u, v) that is a sinc-like spread function. Figs 6A, 6B, and 6C show a selection of possible windows: the standard square window 610 of Fig. 6A, a chequered circular window 620 of Fig. 6B, and a diamond window 630 of Fig. 6C. Other possible windows include, for example, the Bayer pattern 10 window, which models the single colour channels of a digital camera image sensor. In each of the windows 610, 620, and 630, white represents a value I and black a value 0. In Fig. 6A, a square window with unit value 610 is shown. However, there is no restriction to real positive constant window values. In general, complex variable values (with positive or negative components) can be used. The only constraint is that the window transform is function W(u,v) should not be too wildly oscillating, and slowly decaying. The Fourier transform F'(u, v) of the windowed function can now be estimated as follows (including a noise term): F'(u, v)= aoW(u, v)+ ({W(u - u',, v - v,)+ W(u + u', v + v')} 2 + s U -u',v- V',)- W(u ... (19) 2i + N(u,v) 20 A Fourier transform of the padded, window yields the spread function: w(x, y) e "T ,> W(u, v) ... (20) The next step is to estimate the overlap of all the spectral peaks. The manifold interactions are encapsulated by the following matrix equation, which gives the (known) spectral peak values, F', in terms of some (unknown) underlying (sparse) delta 25 spectrum B: 1381845 .DOC IRN: 862327 -24 rF' W, W 2 W1 kK B F' W 2

W

22 W2 .. W2K B 2 _ Wk WK F'= =- l -- W.k K = W.B ...(21) F|' W,,I W WI. W,, . WI B W . f WKI WK 2 WK. WKk WK.. WKK BK The known window interaction of a measurement at position / from a spread function located at k is Wk = W(u, -uk,,v - v) and the complete system is summarised by the matrix W. The known measurements matrix is F'. 5 The solution to the problem, that is to say the underlying pure spectrum, is found by matrix inversion: B = W~'.F' ... (22) For a system with twenty non-zero frequency spectra (N = 20 => K = 2N + 1 = 41 ), a 41 x 41 square matrix inversion is implied. Note that for typical window functions and 10 sparse spectral peaks the matrix is almost diagonal. Closely spaced spectra and/or slowly decaying window transforms cause the matrix inversion to become less stable and harder to invert in the presence of significant noise. In practice, the numerically unstable inversion can be avoided by more appropriate matrix techniques such as Cholesky decomposition, Gauss-Jordan elimination, or Singular Value Decomposition. However, for many 15 applications the above described process can reduce typical spectral estimation error from about 10% to less than 1%. The gain is purely from computational improvements and the data collection (physical measurement) process can be left unaltered. Another significant advantage to the process outlined is that it is computationally efficient. This contrasts dramatically with other approaches based on exhaustive, or iterative, search processes. 20 INDUSTRIAL APPLICABILITY It is apparent from the above that the arrangements described are applicable to the computer and image processing industries. 1381845 1.D)C IRN: 862327 -25 The foregoing describes only some embodiments of the present invention, and modifications and/or changes can be made thereto without departing from the scope and spirit of the invention, the embodiments being illustrative and not restrictive. In the context of this specification, the word "comprising" means "including 5 principally but not necessarily solely" or "having" or "including", and not "consisting only of'. Variations of the word "comprising", such as "comprise" and "comprises" have correspondingly varied meanings. 1381845_1 DOC IRN: 862327

Claims (14)

1. A method of spectral estimation utilising a recorded image of a known test pattern, wherein said test pattern contains a known constellation of spectral peaks, said method comprising the steps of: 5 determining probable peaks in said recorded image; adding said probable peaks to a list of provisional peaks; performing an affine invariant search of all permutations of provisional peaks in said list of provisional peaks using a least squares measure of best fit to identify a best fit; deriving affine parameters from said identified best fit; 10 applying said affine parameters to an ideal constellation to generate a best fit constellation; determining using matrix algebra interaction of all the peaks using a Fourier transform of a window function; disentangling peak values from said determined interactions; and is utilising the disentangled peak values to represent magnitudes and phases of a spectral estimate.

2. The method according to claim 1, wherein a residual between a calculated pattern and a measured pattern indicates a measure of confidence. 20

3. The method according to claim 2, wherein the residual normalised by total energies of both the test pattern image and the recorded image indicates a measure of the confidence. 1381845 1.DOC IRN: 862327 -27

4. The method according to claim 1, wherein the peak detection is implemented as spiral phase detection on a complex derivative of one of the spectral magnitude and intensity. 5

5. The method according to any one of claims I to 4, wherein said constellation defines the number and location of said spectral peaks.

6. The method according to any one of claims I to 5, wherein said recorded image was 1o recorded by an imaging system.

7. The method according to claim 6, wherein the imaging system is chosen from the group of imaging systems consisting of: a digital camera, a document scanner, and a lens of a digital SLR camera. 15

8. The method according to claim 6, wherein said spectral estimate characterises at least one property of said imaging system.

9. A spectral estimation system for determining a spectral estimation of an imaging 20 system by utilising a recorded image of a known test pattern, wherein said recorded image was recorded by said imaging system, said spectral estimation system comprising: an input for receiving said recorded image; a storage medium for storing a provisional list of provisional peaks; a processor for: 25 determining probable peaks in said recorded image; adding said probable peaks to a list of provisional peaks; 13818451 DOC IRN: 862327 -28 performing an affine invariant search of all permutations of provisional peaks in said list of provisional peaks using a least squares measure of best fit to identify a best fit; deriving affine parameters from said identified best fit; s applying said affine parameters to an ideal constellation to generate a best fit constellation; determining using matrix algebra interaction of all the peaks using a Fourier transform of a window function; disentangling peak values from said determined interactions; and 10 an output for presenting said disentangled peak values as representations of magnitudes and phases of a spectral estimate.

10. The spectral estimation system according to claim 10, wherein said imaging system is chosen from the group of imaging systems consisting of: a digital camera, a document is scanner, and a lens of a digital SLR camera.

11. A computer readable storage medium having a computer program recorded thereon, the program being executable by a computer apparatus to make the computer perform a method of spectral estimation utilising a recorded image of a known test pattern, wherein 20 said test pattern contains a known constellation of spectral peaks, said program comprising: code for determining probable peaks in said recorded image; code for adding said probable peaks to a list of provisional peaks; code for performing an affine invariant search of all permutations of provisional peaks in said list of provisional peaks using a least squares measure of best fit to identify a 25 best fit; code for deriving affine parameters from said identified best fit; 1381845_1 DOC IRN: 862327 - 29 code for applying said affine parameters to an ideal constellation to generate a best fit constellation; code for determining using matrix algebra interaction of all the peaks using a Fourier transform of a window function; s code for disentangling peak values from said determined interactions; and code for utilising the disentangled peak values to represent magnitudes and phases of a spectral estimate.

12. A method of spectral estimation utilising a recorded image of a known test pattern, 10 wherein said test pattern contains a known constellation of spectral peaks, said method being substantially as described herein with reference to the accompanying drawings.

13. A spectral estimation system substantially as described herein with reference to the accompanying drawings. 15

14. A computer readable storage medium having a computer program recorded thereon, the program being executable by a computer apparatus to make the computer perform a method of spectral estimation utilising a recorded image of a known test pattern, wherein said test pattern contains a known constellation of spectral peaks, said program being 20 substantially as described herein with reference to the accompanying drawings. DATED this Fourth Day of September, 2008 Canon Kabushiki Kaisha Patent Attorneys for the Applicant 25 SPRUSON & FERGUSON 1381845 .DOC IRN: 862327