Processing methods for coded aperture imaging

Abstract

The present invention relates to processing of coded aperture images. Multiple frames of data acquired by the coded aperture imaging system, each with a different coded aperture array, are processed to form an image. The processing incorporates the constraints that the image solution must be positive and that the solution should be zero outside an expected image region. In one embodiment image enhancement may involve dividing the processed image into image regions having a spatially invariant point spread function and solving an inverse problem for each image region to reduce image blurring.

Description

This invention relates to methods for processing image data acquired in a coded imaging system.

Coded aperture imaging is a known imaging technique which is primarily used in high energy imaging such as X-ray or γ-ray imaging where suitable lens materials do not generally exist, see for instance E. Fenimore and T. M. Cannon, “Coded aperture imaging with uniformly redundant arrays”, Applied Optics, Vol. 17, No. 3, pages 337-347, 1 Feb. 1978. It has also been proposed for three dimensional imaging, see for instance “Tomographical imaging using uniformly redundant arrays” Cannon T M, Fenimore E E, Applied Optics 18, no. 7, p. 1052-1057 (1979)

Coded aperture imaging exploits the same principles as a pinhole camera but instead of having a single small aperture uses a coded aperture mask having an array of apertures. The small size of the apertures results in a high angular resolution and increasing the number of apertures increases the radiation arriving at the detector thus increasing the signal to noise ratio. Each aperture passes an image of the scene to the detector array and so the pattern at the detector array is an overlapping series of images and is not recognisable as the scene. Processing is needed to reconstruct the original scene image from the recorded data.

The reconstruction process requires knowledge of the aperture array used and system configuration and the aperture array chosen is often coded to allow subsequent good quality image reconstruction. The processing is performed using a mathematical model of the particular array at a set location.

Busboom et al. “Coded aperture imaging with multiple measurements” J. Opt. Soc. Am. A, Vol. 14, No. 5, May 1997 propose a coded aperture imaging technique which takes multiple measurements of the scene, each acquired with a different coded aperture array. Busboom et al. discuss that image reconstruction can be performed using a cross correlation technique and, considering quantum noise of the source, discuss the choice of arrays that maximise the signal to noise ratio.

Conventional coded aperture imaging can be thought of as a geometric imaging technique and for the applications it is usually used for, e.g. astronomy, diffraction is negligible. The technique proposed by Busboom et al. discussed above is not concerned with diffraction and hence the effects of diffraction on image reconstruction are not discussed.

Recently, see our co-pending PCT patent application WO2006/125975, it has been proposed to use a reconfigurable coded aperture imager having a reconfigurable coded aperture mask means. The use of a reconfigurable coded aperture mask means allows different coded aperture masks to be displayed at different times. This allows, for example, the direction and FOV of the imaging system to be altered without requiring any moving parts. Further the resolution of the imaging system can also be altered by changing the coded aperture mask displayed on the coded aperture mask means.

The pattern displayed on the coded aperture mask means is a coded aperture mask and at least part of the coded aperture mask is a coded aperture array. That is either the whole pattern displayed on the mask means is a coded aperture array or only part of the pattern is a coded aperture array. For the avoidance of doubt the term aperture used herein does not imply a physical hole in the mask means but merely an area of the pattern which allows radiation to reach the detector.

WO2006/125975 teaches a versatile and lightweight imaging system that can be rapidly configured to have different fields of view or resolution without any moving parts. It eliminates the need for conventional optics, gives conformal imaging capability, can have an infinite depth of field and gives inherent power free encryption since decoding of the image requires knowledge of the coded aperture array used. The imaging apparatus described therein is particularly suitably for several imaging and surveillance applications in the visible, infrared or ultraviolet wavebands.

However, high resolution imaging requires small aperture sizes and a longer optical path from the detector to the mask, which increases the effects of diffraction. Diffraction causes a blurring of the pattern formed by the mask on the detector array, reducing its modulation and having a corresponding detrimental effect on decoded image quality.

It is therefore an object of the invention to provide a method of processing image data from a coded aperture imaging system which recovers good image quality in the presence of diffraction.

Thus according to the present invention there is provided a method of forming an image from data captured from a coded aperture imaging system comprising the steps of taking a plurality of frames of data recorded by the coded aperture imaging system, each frame of data acquired using a different coded aperture array, taking for each frame of data a decoding pattern corresponding to the coded aperture array used to acquire that frame of data and processing the data to provide a decoded image which is a combined regularized least squares fit to all the frames of data. Preferably the processing incorporates at least one of the constraints that the decoded image must be positive and that the decode imaged should be zero outside of an expected image region.

The method according to the present invention therefore uses data acquired using a plurality of different coded aperture arrays, i.e. a plurality of different mask patterns, and attempts to solve the combined inverse problem for the totality of frames of data and corresponding coded aperture arrays.

The skilled person will appreciate that, in the diffraction free case, the signal recorded at the detector array of a coded aperture imaging system can, for a single frame of data recorded with a particular coded aperture array, be described as a convolution of the scene intensity with the mask pattern, i.e. the decoding pattern, of the coded aperture array plus some noise. For the diffraction free case the decoding pattern corresponds to the pattern of the coded aperture array. The object of all decoding algorithms is therefore to recover the scene image by using knowledge of the mask pattern, for instance by performing a deconvolution or cross correlation.

Where diffraction effects are significant however the intensity pattern at the detector array no longer corresponds directly to the aperture function. Instead the diffraction pattern formed at the detector array is in effect a blurred version of the mask pattern. Diffraction effects prevent an exact deconvolution.

The present inventors have found however that using a plurality of frames of data, each acquired using a different coded aperture array, can significantly increase the quality of the recovered image in the presence of diffraction. The method of the present invention may therefore operate on frames of data which are the recorded intensity pattern, formed on the detector of the coded aperture imagine system, for which diffraction effects are significant. By significant is meant that diffraction causes noticeable spreading of the signal and reduces the modulation of the aperture mask.

In the method of the present invention a plurality of sets of data are obtained, each acquired using a different coded aperture array. The method of the present invention seeks to provide a solution which simultaneously decodes the data for each frame data using its associated decoding pattern. As will be described in more detail later, this can be posed as a least squares problem and a solution to the minimisation problem found. It should noted that the method of present invention does not simply combine data from different frames but instead provides a simultaneous solution to the separate problems represented by each frame of data. The decoding pattern used for each coded aperture array preferably corresponds to the diffraction pattern which is formed by the coded aperture array at the detector at the wavelength(s) of interest. As will be described below the diffraction pattern can be measured at the wavelength(s) of operation or can be calculated based on the aperture pattern.

Diffraction effects mean that the least squares minimisation problem is an ill posed problem and therefore some form of regularization is preferably applied to the least squares fit to produce a sensible solution. One preferred form of regularization is a multi-image form of Tikhonov regularization such as proposed by Bertero and Boccacci: M. Bertero and P. Boccacci—“Image Restoration Methods for the Large Binocular Telescope”, Astron. Astrophys. Suppl. Ser. 147, pp. 323-333, 2000. Tikhonov regularisation is a known technique in the solution of inverse problems, see also for example page 108 of “Introduction to Inverse Problems in Imaging, M. Bertero & P. Boccacci, Institute of Physics Publishing, 1998, ISBN 0750304359 (hereinafter referred to as Bertero & Boccacci).

Alternatively a Wiener filter technique could be applied. An iterative technique such as the Landweber iteration could also be used, see page 291 of Bertero & Boccacci.

In order to improve image quality and attempt to overcome diffraction effects the method of the present invention incorporates as much prior information about the solution as possible. Thus the method of the present invention incorporates one or both of the constraints that the image should be positive and that the decoded image should be zero outside an expected region.

Coded aperture imaging is an incoherent imaging technique and the image one is trying to reconstruct is an intensity image. Therefore the true image will not be negative and one can constrain all possible solutions to contain no negative parts. Positively may be incorporated into the solution, for example, using the method described by G D de Villiers, E R Pike and B McNally—“Positive solutions to linear inverse problems”, Inverse Problems, 15, 1999 pp. 615-635.

Conveniently however positivity may be incorporated into the solution using a variant of the Landweber iteration as this is potentially easier to implement (see Bertero and Boccacci, page 291). In this case the method may involve using a projection operator onto the set of functions which is one for all parts of the solution which are positive and zero for all parts of the solution which are negative. In other words only the positive parts of the solution to one iteration are used in the next iteration. Note that the Richardson-Lucy method (also known as the Expectation Maximization method) has similar performance to the projected Landweber method though it is computationally more intensive (Bertero and Boccacci, page 179).

As mentioned above, in the absence of diffraction, coded aperture imaging can be thought of as a geometric imaging technique. Therefore, given a knowledge of the size and spacing of the coded aperture array in the coded aperture imaging system, one knows in advance the expected size and location of the image after it has been decoded. Diffraction effectively causes spreading of the data outside the expected region. The true image should be entirely located within the expected region and hence, ideally, there would be no image outside of this expected region. Therefore a further constraint is that the solution, i.e. the decoded image, must be zero outside an expected region, or in other words the decoded image must have a certain finite support.

The support for the solution, i.e. the constraint that the decoded image must be zero outside an expected region, may also be incorporated into the method of the present invention by using a variant of the Landweber iteration by using a projection operator on the set of functions which is one for all areas within the expected region and zero for all areas outside the region.

If one assumes that the degree of blurring is constant across the image the evaluation step of the Landweber iteration may be performed in the frequency domain as this speeds processing. However to enforce positivity and support for the image one needs to work in the spatial domain. The skilled person will be well aware that one can apply a Fourier transform to the constrained solution to switch to the frequency domain, solve for the new iteration and then apply an inverse Fourier transform to return to the spatial domain.

In an especially preferred method of the present invention therefore the step of processing the data to provide an image solution which is least squares fit to the data from each frame of data comprises an initial step of forming a first estimate of the image solution, and then an iterative improvement step incorporating the constraints that the solution, i.e. the decoded image must be positive and that the solution must be bounded in an expected image region, i.e. that the decoded image should be zero outside of the expected region. The iterative step preferably applies the Landweber interation.

The method of the present invention therefore provides a method of reconstructing good quality images using several different frames of data even when diffraction effects are significant. The method may therefore be used to form images from data acquired in the visible, ultraviolet or infrared spectra where diffraction effects may be significant.

The method may be applied real time using the current frame of data and at least one previous frame of data or it may be applied off line to data previously acquired.

The number of frames of data which are used in the least squares fit may depend on the particular application and in particular on how quickly frames are captured and the evolution of the scene. The method of the present invention is based on the several frames of data being of substantially the same scene. If the scene is changing rapidly or the data is acquired when the coded aperture system is moving rapidly with respect to the scene there may be a limit to how many frames can be combined and still relate to essentially the same scene.

For a static scene combining ten frames or so of data has been found to give good improvement in image quality compared to processing a single frame of data only. The more frames of data corresponding to the same scene that one uses will generally improve image quality but will also increase the processing overhead. Processing m+1 frames will give improved image quality over using m frames but with decreasingly significance as m increases. However the processing load will continue to increase with m and so there will be an optimal number of frames depending on the application. The skilled person will be able to determine the optimal number of frames for a particular application which may be approximately 5 or 10 or 15 or 20 or 25 frames of data or more.

The method of the present invention, as described above, therefore provides benefits compared with processing a single frame of data in producing clearer images and helps overcome diffraction effects. In an embodiment of the invention the data may be further processed to provide further image enhancement. The multiple frame approach described above may therefore form part of an overall method of processing data acquired from a coded aperture imaging apparatus. Such an overall method may take a plurality of frames of data and process each to form an initial image before subsequently performing at least one image enhancement step using the plurality of initial images. The data is processed to form an image for each frame prior to applying the least squares fit.

Each blurred initial image can be seen as the true image convolved with a point spread function. The step of image enhancement is then to recover the true image.

As the skilled person will appreciate the point spread function may vary across the image. The method may therefore involve the step of dividing the each initial image into a series of image regions over which the point spread function is relatively invariant prior to performing the lest squares fit on these image regions to improve quality. Dividing the image into a plurality of small image areas not only ensures that the point spread function is spatially invariant for that area but it also eases the computation as compared with attempting to process the entire image. In processing each small image region the data from each initial image corresponding to that particular region is processed together. Thus the method as described above, using a plurality of different frames of data, is performed separately for each small image region.

The initial image can be formed for each frame of data in a variety of ways but preferably a simple cross correlation is performed for the single frame of data using a decoding pattern appropriate to the coded aperture array used to capture the data. This is the conventional approach to image decoding in coded aperture imaging except that in the presence of diffraction one uses the diffraction pattern of the coded aperture array.

The initial image is therefore decoded using a particular decoding pattern for the whole image. As mentioned above, the decoding pattern corresponds to the diffraction pattern caused by the mask at the wavelength of operation. The skilled person will appreciate however that a point source located at one part of the scene would lead to a different diffraction pattern being formed compared to the same point source located at a different part of the scene. In other words the decoding pattern used in forming the initial image may only be the correct function for part of the scene and hence other parts of the initial image may be blurred due to use of a decoding pattern which is slightly different to the actual decoding pattern that corresponds to that part of the scene. The least squares fit for each image region preferably therefore uses a decoding pattern optimised for that region of the image.

The decoding pattern appropriate to each region could be recorded by imaging a point source having the appropriate spectral characteristics in different regions of the scene and recording the intensity pattern on the detector array for each part of the scene. Our co-pending patent application GB0602373 teaches various methods by which the appropriate decoding pattern could be recorded. Alternatively the decoding pattern for each image region could be calculated.

Preferably the step of processing each image region, the least squares fit, involves a Tikhonov regularization. Preferably the step of processing each image region also incorporates the constraint that the solution should be positive.

Processing discrete image regions could lead to image artifacts at the edge of adjacent regions. Therefore preferably a first image region is taken and processed. However only the values of the reconstruction, i.e. decoded image formed by the least squares fit, at the centre of this region are stored. The region is then moved by a few pixels (so that the second image region overlaps with the first image region) and the process is repeated. This avoids any image artifacts being created.

The present invention therefore can be implemented as a two stage process. Thus in another aspect of the invention there is provided a method of forming an image from data captured from a coded aperture imaging array comprising the steps of taking a plurality of frames of data, each frame of data acquired using a different coded aperture array, taking for each frame of data a decoding pattern corresponding to the appropriate coded aperture array and processing the data to provide a first image solution for each frame of data and then deblurring the plurality of first image solutions using a set of deblurring functions derived from a set of decoding patterns each optimised for different regions of the image. Preferably the step of deblurring the first image solution incorporates the constraint that the solution must be positive.

It should be noted that whilst the present invention relates to use of multiple mask images the step of dividing a first image solution into a series of sub-images and processing each separately can be advantageously applied even with only a single frame of data. In this case a single frame of data could be processed to form a first image solution using a particular decoding pattern and then divided into a plurality of image regions, each image region being processed separately using a decoding pattern optimised for that image region. This processing of the separate image regions with optimised decoding patterns improves the quality of the image as compared with the first image solution.

It would therefore be possible to perform a method of forming an image from data recorded by a coded aperture imaging system wherein the method comprises the step of dividing the image into a plurality of smaller image regions and processing each small image region separately using a decoding pattern optimised for that image region. In effect one could solve the inverse problem separately for each image region.

Solving the inverse problem for each image region can proceed using the method outlined above, i.e. the initial solution can be solved with a Tikhonov regularisation using Fourier methods. Alternatively a truncated singular function expansion could be used as proposed in “Scanning singular-value decomposition method for restoration of images with space-variant blur”, D A Fish, J Grochmalicki and E R Pike, J. Opt. Soc. Am A, 13, no. 3 1996, pp 464-469. This is a more computationally intense method that Tikhonov regularisation. This method requires calculation of the singular value decomposition (SVD) associated with the point-spread function. However if the point-spread function is spatially invariant then the SVD just needs to be calculated once for treating the entire image. Again one would preferably process a plurality of overlapping small image regions to form a solution for each region and would retain the values for the particular solution at the centre of each region only so as to avoid any edge artifacts.

In general other prior information about the scene or targets in the scene can be incorporated into the method of processing the data to form an image. For instance, where one has prior knowledge that the image has a small number of point targets the method may involve an image enhancement step may that uses super-resolution methods involving curve-fitting.

Further, if the coded aperture imager is to be used for tracking then high-resolution patches may only be needed where tracking is being carried out. This would cut the computational load significantly. Therefore the method may perform image enhancement only at parts of interest of the image, i.e. moving parts of the scene or areas of possible or confirmed targets.

Where the coded aperture imager is used for target tracking information from more than one image can also be combined in a track-before-detect scheme. Track Before Detect algorithms have been previously used in the radar and sonar fields and use data from several acquisitions of data from the sensor together to improve target identification. A similar approach could be used to different images from a coded aperture imaging system.

The foregoing description has described the invention in terms of a method only. However the present invention also allows a novel coded aperture imaging system. Hence, in another aspect of the invention there is provided a coded aperture imaging system comprising a detector array arranged to image a scene through a reconfigurable coded aperture mask means and a processor, wherein the reconfigurable coded aperture mask means is adapted to provide a plurality of different coded aperture arrays in series and the processor records the signal from the detector array for each different coded aperture array and applies the method as described above to form an image.

The present invention will now be described by way of example only with respect to the following drawing of which;

FIG. 1 shows schematically a coded imaging system of the present invention,

FIG. 2 illustrates how a decoding pattern for each coded aperture can be recorded,

FIG. 3 illustrates recording a decoding pattern for each coded aperture array for a different part of the scene,

FIG. 4 illustrates acquiring a plurality of frames of data using a coded aperture imager,

FIG. 5 illustrates the a first stage of processing of the data according to an embodiment of the present invention,

FIG. 6 illustrates a second stage of processing of the data in an embodiment of the invention

FIG. 7 illustrates in more detail the step of forming a plurality of local deblurring functions as shown in FIG. 6,

FIG. 8 illustrates in more detail the step of local multimask processing shown in FIG. 6,

FIG. 9 illustrates how the method of the present invention operates a) in the case of processing a whole image and b) in processing separate image regions, and

FIG. 10 shows two example processed images.

Coded aperture imaging (CAI) is based on the same principles as a pinhole camera. In a pinhole camera, images free from chromatic aberration are formed at all distances away from the pinhole, allowing the prospect of more compact imaging systems, with a much larger depth of field. However, the major penalty is the poor intensity throughput, which results from the small light gathering characteristics of the pinhole. Nevertheless, the camera is still able to produce images with a resolution determined by the diameter of the pinhole, although diffraction effects have to be considered. The light throughput of the system can be increased by several orders of magnitude, while preserving angular resolution, by using an array of pinholes. Each detector element sees the result of the summation of contributions from the various pinholes, corresponding to each viewpoint of the scene.

Another way of understanding the operating principle of conventional CAI is to observe that this is a purely geometric imaging technique. Light from every point in a scene within the field of regard (FOR) of the system casts a shadow of the coded aperture onto the detector array. The detector measures the intensity sum of these shadows. The coded aperture is specially designed such that its autocorrelation function is sharp with very low sidelobes. Typically pseudorandom or uniformly redundant arrays (URA) (such as described in E. Fenimore and T. M. Cannon, “Coded aperture imaging with uniformly redundant arrays”, Applied Optics, Vol. 17, No. 3, pages 337-347, 1 Feb. 1978) are used where a deconvolution or decorrelation of the detector intensity pattern can yield a good approximation to the point distribution in the scene.

Previous CAI systems have generally been employed in application where diffraction effects are minimal. For instance coded aperture imaging has often been used in astronomical imaging. However, for some applications of coded aperture imaging techniques, it is necessary to improve the angular resolution significantly. This can be especially true when operating in the visible, infrared or ultraviolet wavebands say, or in other wavebands requiring high resolution imagery. Assuming that the detector pixels are smaller than the feature size of the coded aperture array, p, the angular resolution is determined by tan⁻¹ (p/s) where s is the optical distance between the mask and the detector array. Therefore increasing the resolution of the imager requires either decreasing the size of the apertures or increasing the mask to detector distance or both. With relatively small apertures and/or large mask to detector distances, diffractive effects start to become significant. The blurring effects of diffraction mean that the pattern projected onto the detector array is effectively smeared, reducing image quality.

FIG. 1 shows schematically an example of a coded aperture imaging system, generally indicated 2. Rays of light from points in the scene 4 fall onto a particular coded aperture array 6. The coded aperture array acts as a shadow mask and therefore a series of overlapping coded images are produced on the detector array 8. At each pixel on the detector array, the intensities from the overlapping, coded images are summed. The output from the detector array 8 is passed to a processor 10 where image of the scene can be subsequently decoded from the detector signals using a variety of digital signal processing techniques.

Recently WO2006/125975 has proposed using a reconfigurable mask means 6 to provide a reconfigurable coded aperture array. The coded aperture mask means 6 is controlled by a controller 12 which controls the reconfigurable mask means to display different coded aperture masks. If only part of the coded aperture mask means displays a coded aperture array, the rest of the mask preventing radiation from reaching the detector, then the field of view of the apparatus is determined by the location and size of the coded aperture array relative to the detector, changing its position or size changes the field of view and/or resolution of the imager.

CAI therefore offers the ability to provide a compact and lightweight imaging apparatus which has a large depth of field and has a changeable field of view without requiring any moving parts.

As mentioned above the intensity pattern formed at the detector array 8 is a series of overlapping images of the scene created by the various apertures. This intensity pattern requires decoding to construct a recognisable image of the scene.

FIG. 4 illustrates the capture of multiple frames of data. Reconfigurable coded aperture mask means 6 is sequentially configured to provide a series of coded aperture masks 20 and the output 22 of the detector array 8 is recorded for each coded aperture mask displayed.

For a single mask pattern the basic equation relating the recorded data to the desired scene isg=Af+noise  eqn (1)where g is the recorded data, A is a linear operator containing the diffraction pattern of the mask and f is the desired scene. The camera will only record intensities so, for the case of narrowband radiation, A will represent the modulus squared of the diffraction pattern of the mask, as recorded on the detector array, and the operator A will be a convolution operator. At this stage it is assumed that the detected radiation is narrowband, for instance, as would be the case using a narrowband detector or applying a narrow band filter.

The multi-mask approach involves decoding simultaneously sets of data corresponding to different mask configurations, e.g. M₁-M₄. We are required to simultaneously solve for f a set of operator equations of the form:A_if=g_j, j=1, . . . , p  eqn (2)where p is the number of different mask configurations used. This can be phrased as a least-squares problem, as described in M. Bertero and P. Boccacci—“Image Restoration Methods for the Large Binocular Telescope”, Astron. Astrophys. Suppl. Ser. 147, pp. 323-333, 2000. where one seeks to minimise over f the discrepancy functional

$\begin{array}{cc}{\varepsilon }^2\left(f\right)=\sum _{j=1}^p{A_{j}f-g_j}^2& \mathrm{eqn}\left(3\right)\end{array}$

The solution to this minimisation problem, f, will satisfy the normal equation:

$\begin{array}{cc}\sum _{j=1}^p{A}_j^T{A}_{j}f=\sum _{j=1}^p{A}_j^T{g}_j& \mathrm{eqn}\left(4\right)\end{array}$

Now noting that the operators A_j (at least for the narrowband problem) are convolution operators we can take the two-dimensional DFT of both sides of this equation to yield

$\begin{array}{cc}\sum _{j=1}^p{{\hat{A}}_j\left(\omega \right)}^2\hat{f}\left(\omega \right)=\sum _{j=1}^p{\hat{A}}_j^T\left(\omega \right){\hat{g}}_j\left(\omega \right)& \mathrm{eqn}\left(5\right)\end{array}$ where Â_j^T(ω) is the conjugate of the Fourier transform of the diffraction pattern of the jth mask and ĝ_j(ω) is the Fourier transform of the data corresponding to the jth mask. This is the key equation one needs to solve.

In the presence of diffraction solution of the basic equation is an ill-posed inverse problem and some form of regularisation is necessary to produce a sensible solution. One possibility is a multi-image form of Tikhonov regularisation (Bertero and Boccacci, 2000, equation 13).

In the Fourier domain the Tikhonov solution takes the form

$\begin{array}{cc}\hat{f}\left(\omega \right)=\sum _{j=1}^p\frac{{\hat{A}}_j^T\left(\omega \right)}{{\hat{A}\left(\omega \right)}^2+\mu }{\hat{g}}_j\left(\omega \right)& \mathrm{eqn}\left(6\right)\end{array}$ where μ is the regularisation parameter and

$\begin{array}{cc}{\hat{A}\left(\omega \right)}^2=\sum _{j=1}^p{{\hat{A}}_j\left(\omega \right)}^2& \mathrm{eqn}\left(7\right)\end{array}$

This solution does not take into account the prior knowledge that the solution must be non-negative and must have a certain, finite, support, i.e. one knows where the image is expected to be in image space. A simple way of doing this is to modify the usual Landweber iteration.

The method described here is a modified version of the Landweber iteration given in eqn (17) in M. Bertero and P. Boccacci—“Image Restoration Methods for the Large Binocular Telescope”, Astron. Astrophys. Suppl. Ser. 147, pp. 323-333, 2000. The modification is necessary in order to incorporate into the problem knowledge of the support of the image.

The iteration takes the form

$\begin{array}{cc}{f}_{+}^{k+1}=P_{+}\left\{f_{+}^k+\tau \sum _{j=1}^p\left(A_j^T{g}_j-A_j^T{A}_j{f}_{+}^k\right)\right\}& \mathrm{eqn}\left(7\right)\end{array}$

Here k is the iteration number and τ is the relaxation parameter. The + subscript on f indicates that the eventual solution will be non-negative. P₊ is the projection operator onto the set of non-negative functions, defined by

$\begin{array}{cc}\left(P_{+}f\right)\left(x\right)=\begin{array}{c}f\left(x\right),\\ 0,\end{array}\begin{array}{c}f\left(x\right)\ge 0\\ f\left(x\right)<0\end{array}& \mathrm{eqn}\left(8\right)\end{array}$

In addition to this projection, the projection onto the set of functions with a given support S, i.e. the expected part of image space where the image will be formed, needs to be carried out. This is accomplished using the projection operator P_S defined by

$\begin{array}{cc}\left(P_{S}f\right)\left(x\right)=\begin{array}{cc}f\left(x\right),& x\in S\\ 0,& x\notin S\end{array}& \mathrm{eqn}\left(9\right)\end{array}$

For the problem in hand, due to the size of the images, it is much easier to carry out the convolutions in the Fourier domain. The iteration would then proceed by swapping between the spatial domain and the Fourier domain, with the projection onto the set of non-negative functions and the projection onto the set of images with a certain support being carried out in the spatial domain.

The method of the present invention therefore generally proceeds as follows. First an the initial estimate of the solution is constructed using multi-image Tikhonov regularisation (equation 6 above). This result is then projected onto the set of images with a given support (i.e. set all elements of the image outside the known region of interest to zero).

It should be noted that the outcome of the solution to the Tikhonov regulation is not constrained to be negative and generally will have significant negative as well as positive parts. One could simply set any negative parts of this solution to zero and progress using the positive parts. However it has been found that applying such a technique and then performing the iteration described above to a greyscale image, as the iteration progressed, the shadows on the image tended to increase in size, giving an unnatural appearance to the image. It was determined that this is due to the taking of the positive part of an initial image which has significant negative parts. The present inventors have found that with this type of starting point once the negative regions on the initial image have been set equal to zero they tend, as the iteration proceeds, to remain zero which can distort the final solution.

Since the initial image for the iteration stage does not need to agree accurately with the data (the agreement being reached through the iteration) one can shift the initial image to make it non-negative. This seems to give very beneficial results. Therefore a data shift is applied to this initial solution to create a solution which is entirely positive and this is stored as the initial f_new, The data shift is in effect a simple addition to each pixel of a value which makes every pixel non-negative.

Having thus established an initial solution f_new the main iteration is repeated a predetermined number of times, the main iteration comprising the steps of:

• (a) Project f_new onto the known support using P_S and store the result in f_new
• (b) Take the positive part of f_new using P₊ and store the result in f_new
• (c) Take the DFT of f_new, Call this {circumflex over (f)}_old
• (d) Evaluate{circumflex over (f)}_new(ω)=τ(A^Tg)(ω)+(1−τ(A^TA)(ω)){circumflex over (f)}_old(ω)  eqn (10)where:

$\begin{array}{cc}\left(A^{T}g\right)\left(\omega \right)=\sum _{j=1}^p{\hat{A}}_j^T\left(\omega \right){\hat{g}}_j\left(\omega \right)\mathrm{and}& \mathrm{eqn}\left(11\right)\\ \left(A^{T}A\right)\left(\omega \right)=\sum _{j=1}^p{{\hat{A}}_j\left(\omega \right)}^2.& \mathrm{eqn}\left(12\right)\end{array}$

• (e) Take the inverse DFT of {circumflex over (f)}_new:f_new=F⁻¹({circumflex over (f)}_new)  eqn (13)

The relaxation parameter τ is chosen to satisfy

$\begin{array}{cc}0<\tau <\frac{2}{{\hat{A}}_{\max }^2}\mathrm{where}& \mathrm{eqn}\left(14\right)\\ {\hat{A}}_{\max }^2=\begin{array}{c}\max \\ \omega \end{array}\sum _{j=1}^p{\hat{A}}_j^T\left(\omega \right)·{\hat{A}}_j\left(\omega \right)& \mathrm{eqn}\left(15\right)\end{array}$

The iteration described above has been used on various data sets and a relaxation parameter of;

$\begin{array}{cc}\tau =\frac{1.8}{{\hat{A}}_{\max }^2}& \mathrm{eqn}\left(16\right)\end{array}$ has been found to give good results.

As mentioned above the method described above has been derived based on the assumption of narrow band radiation. The method can also be used in a broad band situation if the spectrum of radiation from each point in the scene is the same as at any other point. In this case the problem simplifies back to a convolution equation and hence the previous mathematics can still be used. The basic equation in the spatial domain can be written asg(y)=∫∫K_λ(y−x)f(λ,x)dxdλ  eqn (17)where K_λ(x) is the diffraction pattern of the mask at wavelength λ. If one assumes that the spectrum of light coming from any point of the scene is the same as that at any other point then one may writef(λ,x)=h(λ)f(x)  eqn (18)

Interchanging the order of integration in eqn (17) then yieldsg(y)=∫└∫K_λ(y−x)h(λ)dλ┘f(x)dx  eqn (19)

The term in brackets is still a function of (y−x). If one denotes this function by {tilde over (K)}(y−x) then one has:g(y)=∫{tilde over (K)}(y−x)f(x)dx  eqn (20)

This is still a convolution equation and so the methods devised for the narrowband convolution problem can all be used as they stand. The only difference is that one must use the broadband point-spread function {tilde over (K)} instead of K_λ.

FIG. 9 a illustrates the principles of the invention. For each frame of data 22 is taken an appropriate decoding pattern 24 which are then processed, step 26, using the method described above to yield an image. The decoding patterns may be calculated for each mask pattern or, as shown in FIG. 2, could be actually measured by the coded aperture imaging system.

As shown in FIG. 2 a point source 30 may be located in the field of view of the imager. The point source is chosen to match the expected spectral characteristics of the scene so that it forms an appropriate diffraction pattern. The output from the detector array for each coded aperture mask 20 is recorded and represents the decoding pattern 24 for that mask.

The present invention therefore uses a plurality of frames of data to form an image which is clearer than could be formed using a single frame of data alone. In another embodiment of the invention however further clarity can be obtained though by processing to account for local blurring across the image.

In this embodiment a plurality of local deblurring functions, i.e. deblurring functions appropriate for a small image region, are obtained by using a decoding pattern which is optimised for each image region. Again the decoding pattern can be recorded as described above with reference to FIG. 2. FIG. 3 shows a point source located at a different part of the scene. Again the plurality of coded aperture mask patterns 20 are sequentially written to the reconfigurable coded aperture array means and the output of detector array 8 recorded for each coded aperture array. This provides a series of decoding patterns 34 which are optimised for the part of the image which corresponds to the location of point source 32 in the scene. This step is repeated by moving the point source relative to the coded aperture imager so that it appears at a plurality of different positions in the scene, n different positions, and recording a decoding pattern for each coded aperture array for each position.

Having recorded (or calculated) the decoding pattern for each coded aperture array for each position in the scene the coded aperture imager is arranged to view a scene as described above and shown in FIG. 4 and record a plurality of frames of data 22.

In this embodiment however each frame of data 22 is initially correlated, see FIG. 5, with an appropriate decoding pattern 24 to yield a partially deblurred image 28 for each frame of data. It is convenient to use the set of decoding patterns 24 which correspond to the central image region, i.e. those that would be recorded using a point source which is directly in front of the coded aperture imager (which will be referred to as an on axis point source and the corresponding decoding pattern as an on axis decoding pattern), however the set of decoding patterns corresponding to any one image region could be used. This results in a set of m initial images, each of which is at least partially deblurred but for which the point spread function will vary over the image. The next stage of processing therefore divides each of the m images into n separate image regions, over which the point spread function is spatially invariant, and deblurs each image region separately.

Referring now to FIGS. 6 and 7 the whole set of decoding patterns 24, 34 is taken and used to form a set 36 of n local deblurring functions or point spread functions (see FIG. 7), one for each image region. The point spread functions to be deconvolved in this further deblurring step are calculated, as shown in FIG. 7, by forming the cross-correlation of the on-axis decoding pattern (or whatever set of decoding patterns were used to form the initial image) with the decoding pattern for the appropriate image region, for example position 2 shown in FIG. 3.

In other words the off axis point spread functions are calculated by decoding the mask diffraction pattern obtained when the coded aperture imager is illuminated with an off axis point source. The decoding is performed using a decoding pattern derived from the mask diffraction pattern when the system is illuminated with an on axis point source.

Referring now to FIG. 8, the m partially partially deblurred images are passed to an image segmentation algorithm 80. This divides each image into n different image regions. This therefore results in a set of m by n partially deblurred image regions 82.

For each region, in general the ith region 84, the appropriate part of each of the m partially deblurred images is taken along with the corresponding set of local deblurring functions 86.

Thereafter the method as described above can be applied to each image region. In other words a Tikhonov regularisation 38 is applied to the combined problem and then a Landweber iteration performed incorporatinges the constraint that the image must be positive and should be zero in regions of the image that are known to be dark in the scene. This results in a deblurred image for the ith region. The method is repeated for each image region. The n image regions processed overlap and the results of processing for each image are only retained in the centre of the image region to avoid any edge artifacts. This results in a image 40 which is based on the combined frames of data and which is deblurred over the n different image regions.

As mentioned above the same processing is applied to each image region in the embodiment which uses local deblurring functions as is applied in the first embodiment without local deblurring. There is a slight difference in what the inputs to the method are, as illustrated in FIG. 9, but the principle is the same. FIG. 9a shows the inputs for the first embodiment of the invention where the whole image is processed. In this case the inputs are the m intensity patterns as recorded by the detector array and the appropriate m decoding patterns. These are used to form a combined inverse problem which is solved by a regularised least squares fit.

In the second embodiment shown in FIG. 9b, where local deblurring functions are also used, each image region is processed separately. In order to be able to select each image region each frame of data is initially processed to form a first image using, usually, the on axis decoding pattern. The inputs to the multi-frame processing therefore comprise a plurality of sub-images, i.e. the same part of each of the m first images, and also an appropriate deblurring function for each sub-image. The problem to be solved is a deconvolution, exactly as before in the first embodiment and so the skilled person will understand that the same processing method can be applied.

FIGS. 10 a and 10b show the improvement due to multimask processing in two processed images of the same test scene which is a standard resolution chart. FIG. 10a shows the processed image using a single mask and Landweber iteration. FIG. 10b shows the final image from multimask processing using eight masks and multimask Landweber iteration. It can be seen that the multimask processing results in a much clearer final image.

Claims

1. A method of forming an image from data captured from a coded aperture imaging system comprising the steps of taking a plurality of frames of data recorded by the coded aperture imaging system, each frame of data acquired using a different coded aperture array, taking for each frame of data a decoding pattern corresponding to the coded aperture array used to acquire that frame of data and processing the data to provide an image solution which is a combined regularised least squares fit to all the frames of data.

2. A method of forming an image as claimed in claim 1 wherein the processing incorporates at least one of the constraints that the solution must be positive and that the solution must be zero outside of an expected image region.

3. A method of forming an image as claimed in claim 1 wherein the data recorded by the coded aperture system is an intensity pattern of radiation for which diffraction effects are significant.

4. A method as claimed in claim 1 comprising a multi-image Tikhonov regularization.

5. A method as claimed in claim 1 wherein the step of processing the data incorporates the constraint that the solution is positive using a Landweber iteration.

6. A method as claimed in claim 5 wherein the Landweber iteration applies a projection operator onto the set of functions which is one for all parts of the solution which are positive and zero for all parts of the solution which are negative.

7. A method as claimed in claim 1 wherein the step of processing the data incorporates the constraint that the solution is zero outside an expected image region using a Landweber iteration.

8. A method as claimed in claim 6 wherein the Landweber iteration applies a projection operator onto the set of functions which is one for all areas within the expected image region and zero for all areas outside the expected image region.

9. A method as claimed in claim 5 wherein the evaluation step of the Landweber iteration is performed in the frequency domain and the constraints on the solution are imposed in the spatial domain.

10. A method as claimed in claim 1 wherein the step of processing the data to provide an image solution which is a least squares fit to the data comprises an initial step of forming a first estimate of the image solution, and then an iterative improvement step incorporating the constraints that the solution must be positive and that the solution must be zero outside an expected image region.

11. A method as claimed in claim 10 wherein the iterative step applies the Landweber iteration.

12. A method as claimed in claim 1 comprising an initial step of processing each frame of data to form an initial image prior to performing the least squares fit.

13. A method as claimed in claim 12 further comprising the step of dividing each initial image into a series of image regions prior to performing the least squares fit.

14. A method as claimed in claim 13 wherein a separate least squares fit is applied to the data for each image region.

15. A method as claimed in claim 14 wherein the least squares fit for each region uses a decoding pattern optimised for that image region.

16. A method as claimed in claim 14 wherein the least squares fit to the data for each image region comprises a Tikhonov regularization.

17. A method as claimed in claim 14 wherein the least squares fit to the data for each image region incorporates the constraint that the solution should be positive.

18. A method as claimed in claim 14 wherein the image regions overlap and only the values of the least squares fit at the centre of each region are stored.

19. A method as claimed in claim 12 wherein the step of processing each frame of data to form an initial image comprises performing a cross correlation on each frame of data with a decoding pattern appropriate to the coded aperture array used to capture that frame of data.

20. A method of forming an image from data captured from a coded aperture imaging array comprising the steps of taking a plurality of frames of data, each frame of data acquired using a different coded aperture array, taking for each frame of data a decoding pattern corresponding to the appropriate coded aperture array and processing the data to provide a first image solution for each frame of data and then deblurring the plurality of first image solutions using a set of deblurring functions derived from a set of decoding patterns each optimised for different regions of the image.

21. A method as claimed in claim 20 wherein the step of deblurring the first image solution incorporates the constraint that the solution is positive.

22. A method as claimed in claim 1 wherein the data recorded by the coded aperture imaging system is the intensity pattern of radiation within the visible, ultraviolet or infrared spectra.

23. A coded aperture imaging system comprising a detector array arranged to image a scene through a reconfigurable coded aperture mask means and a processor, wherein the reconfigurable coded aperture mask means is adapted to provide a plurality of different coded aperture arrays in series and the processor records the signal from the detector array for each different coded aperture array and applies the method as claimed in claim 1 to form an image.