Transcript Slide 1
The Large Binocular Telescope: a laboratory for image reconstruction
M. Bertero DISI – Universita’ di Genova
Pisa
–
October 12, 2005
Joint work with:
- Barbara Anconelli - Genova - Patrizia Boccacci - Genova - Marcel Carbillet - Nice - Serge Correia - Potsdam .Gabriele Desiderà - Genova - Henri Lanteri - Nice
Outline
- The Large Binocular Telescope (LBT) Restoration of LBT images - The standard approach - Correction for boundary effects - Super-resolution A general approach Objects with high dynamic range Concluding remarks
The Large Binocular Telescope (LBT) – 1/7
The LBT, under construction on the top of Mount Graham, Arizona, is a new-conception telescope, consisting of: - two 8.4 m mirrors, with the centers at a distance of 14.4 m; - a very advanced multi-conjugate AO system (adaptive secondaries, pyramid sensors, etc.); One of the instruments will be a Fizeau interferometer for the combination of the images of the two mirrors.
http://lbtwww.arcetri.astro.it
First mirror already installed and alluminized Second mirror transported from the Mirror Lab to Mount Graham First interferometric light : 2007 – 2008 (?)
The Large Binocular Telescope (LBT) - 2/7
The structure pre-assembled in Milan, June 2001 The enclosure on the top of Mount Graham, Arizona
The Large Binocular Telescope (LBT) - 3/7
LINC/NIRVANA (Lbt Interferometric Camera / Near-IR Visible Adaptive iNterfermeter for Atronomy) is the German Italian beam combiner for LBT.
Will combine the radiation from the two 8.4 m primary mirrors in a “Fizeau” mode, allowing true imagery.
Will operate at wavelengths between 0.6 and 2.4 microns.
The first camera will cover 10 arcsec with a pixel size of about 5 mas in K-band (2048 x 2048 pixels).
Multiconjugate AO correction is over a field of about 2 arcmin so that even larger images will be possible in the future (24000 x 24000 ?).
The Large Binocular Telescope (LBT) – 4/7
In case of perfect circular mirrors, perfect interferometer, no atmosphere, the Point Spread Function (PSF) of LBT is given by (we assume that the telescope is aligned in the x -direction):
K
( x , ) 2 1 2
J
1 ( 1 1 | | | |) 2 cos 2 ( 2 x ) , | | x 2 2 ; 1 2
D
,
D
diameter of the mirrors 8.4
m
B
,
B
distance between mirrors centers 14.4
m First zero in the x direction : x 1 = /2B ; First zero in the direction : 1 = 1.22 /D .
The Large Binocular Telescope – 5/7
x, plane u,v - plane
The Large Binocular Telescope – 6/7
8.4 m 22.8 m
Imaging with LINC/NIRVANA : with a few observations at different parallactic angles and suitable reconstruction algorithm, the image, as seen by
a 22.8-m telescope,
can be reconstructed !…
The Large Binocular Telescope – 7/7
* = * = * = incoherent sum NO !
Restoration of LBT images – 1/4
In the case of LINC/NIRVANA the image restoration problem consists in estimating a unique high-resolution image from p images g j , corresponding to p different orientations of the baseline, being given the corresponding PSFs K j and backgrounds b j . Therefore it is a problem of multiple images deconvolution. Accuracy and resolution are strongly related to the coverage in the u,v plane.
In astronomy the most frequently used restoration method is an iterative method known as RL-method (Richardson, 1972; Lucy, 1974); in the case of deconvolution problems, it coincides with the ML-EM method (Shepp and Vardi, 1982), proposed in tomography.
Restoration of LBT images – 2/4
Model for the p images acquired by LINC/NIRVANA CCD camera (based on Snyder et al, 1993): for j=1,....,p g j (n) g j, s (n) g j, b (n) r j (n) g j,s (n) = number of photons coming from the source; Poisson process with expected value (A j f)(n) = (K j *f)(n) (K j =PSF of the j th baseline orientation, f(n’) = expected value of the number of photons emitted at pixel n’); g j,b (n) = number of photons due background emission, etc; Poisson process with (constant) expected value b j ; r j (n) = read-out noise due to the amplifier; realization of an independent Gaussian process with expected value c and standard deviation s .
Restoration of LBT images – 3/4
Summary of the main features of astronomical images - Bandlimiting (the support of the FT in the {u,v}-plane is bounded) Background due to sky emission (to be estimated) - Contamination by both Poisson and Gaussian noise (whose statistical parameters can be estimated) - High Dynamic range (very often orders of magnitudes between the brightness of different objects in the field of view)
Restoration of LBT images – 4/4
Properties to be satisfied by the restored images: Nonnegativity (assuming an estimate of the background) Conservation of the total flux Correct estimation of the relative local fluxes magnitude of stars: photometry) (relative Correct estimation of the relative angular separation (astrometry)
The standard approach – 1/5
In the case of images dominated by photon noise, the read out noise can be neglected and the likelihood function is given by: P(g | f) j p 1 n e [(A j f)(n) b] [(A j f)(n) b j ] g j (n) g j (n) !
To maximize this function is equivalent to minimize the Kullback-Leibler directed distance or Csisz ár I-divergence measure of the discrepancy between the detected images g j and the computed images A j f + b j (j=1,...,p):
J
0
( Af ; g )
j p
1 n
g
j
(n) ln (A
j
g
j
(n) f)(n)
b
j
(A
j
f)(n)
b
j
g
j
(n)
.
The standard approach – 2/5
If we denote by A j T method becomes: the transposed matrices
,
the standard EM - initialize with f (0) >= 0; - given f (k) , compute f (k+1) by:
f ( k
1
)
1
p f ( k) j p
1
A j T A j f g j ( k)
b j
.
Again it is assumed that each PSF is normalized in such a way that the sum of its pixel values is one.
The method is slow and computationally heavy.
The standard approach – 3/5
The iterates converge, for any initial guess, to a minimum point of the I-divergence; the minimum may not be unique and may be corrupted by strong noise propagation.
Each iterate is non-negative.
The total flux is correctly reproduced (exactly in the case of zero background).
In the case of sparse and resolved objects (systems of isolated stars at an angular distance greater than the diffraction limit) the photometric (relative magnitudes) and astrometric (relative angular separations) parameters are correctly given.
In the case of diffuse objects the iterates have a semi convergent behaviour.
The standard approach – 4/5
The analogy between LBT-images and projections in tomography suggests the extension to LBT (Bertero and Boccacci, 2000) of the OSEM (Ordered Subsets Expectation Maximization) method (Hudson and Larkin, 1994): initialize with
f (
0
)
0 , satisfying the flux constraint , given
f ( k)
, set
j
k
1 mod
p
and compute :
f ( k
1
)
f ( k) A j T A j f g j ( k)
b j
; renormaliz e
f ( k
1
)
in such a way that it satisfies the flux constraint .
The computational cost per iteration is reduced by a factor of about 4/(3p+1).
The standard approach – 5/5
Reconstruction of a diffuse object
original object image with PSF @ 60 deg.
reconstruction 0.7”
Correction for boundary effects – 1/4
When the object extends beyond the boundary of the FOV, the FFT-based methods are equivalent to use periodic boundary conditions; hence they introduce discontinuities at the boundary, which, in the deconvolved image, produce Gibbs oscillations, also called
ripples
.
A simple method has been proposed, which consists in reconstructing the object over a domain broader than the FOV, hence without introducing artificial boundary conditions.
One can apply the Richardson-Lucy (or OSEM) method to this problem; the use of FFT is still possible using arrays 2Nx2N and extending the detected images by zero padding.
The results are quite satisfactory.
Correction for boundary effects – 2/4
Scheme of the approach :: bottom left: the array of the RLM a factors bottom right: the windowed weights, defining the reconstruction region
Correction for boundary effects – 3/4
The modified OSEM algorithm is as follows (the bar denotes 2Nx2N arrays: initialize compute with a
j f
(
0
)
A T j
0 , satisfying
M S
(
w R
,
j
1 / a
j j
the flux 1 ,..., ,
n
R
;
w R
,
j
constraint
p
) , 0 , elsewhere; given
f ( k)
, set
j
k
1 mod
p
and compute : -
f ( k
1
)
renormaliz e
w f R
,
j f ( k)
A j T A j f g ( k) j
b j
;
( k
1
)
in such a way that it satisfies the flux constraint .
OSEM
Correction for boundary effects - 4/4
Ideal PSFs SR=70% SR=26% New method
Super-resolution – 1/6
• In astronomy the term denotes any method able to provides a resolution better than the diffraction limit; • if is the wavelenght the limit is: /D in the case of a single mirror with diameter D; /B in the case of LBT with maximum baseline B.
• the amount of super-resolution is controlled by: – the angular Space – Bandwidth product (SBP), given by the ratio between the angular size of the object and the diffraction limit; – Signal-to-noise ratio (SNR): the amount of super resolution increases with increasing SNR
Super-resolution – 2/6
The standard restoration method implements the constraint of non-negativity as well as a constraint on the total flux of the object and therefore it already provides a moderate amount of super-resolution.
However the crucial constraint is on the support of the object. It can be implemented by using a ’’localization’’ property of the algorithm: due to the ’’multiplicative’’ structure of the algorithm, if the initial guess is zero in one pixel, all the iterates are zero in that pixel.
If the domain of the object is known, one can initialize the algorithm with the characteristic function of the domain.
Super-resolution – 3/6
• construction of a mask: matrix with values 0 or 1: – 1 where the image resulting from the first step is greater than a fixed treshold = given percent of the image maximum value – 0 elsewhere.
Result of the first step Mask
Super-resolution – 4/6
2. RL restoration initialized with the mask of the domain as estimated in step 1.
• • smaller number of iterations (250-5000) can already provide a super-resolved image but the photometry may not be correct OS-EM 5000 iter
Result of the first step Result of the second step
Super-resolution – 5/6
If the photometry is not correct the first two steps are used to estimate the positions of the two stars
3.
solution of a least-squares problem
,
where the unknowns are the magnitudes of the stars original difference of magnitude difference of magnitude after step 3
Super-resolution – 6/6
Minimum separation limit in the case of binaries with the same magnitude number of iterations: Step 1: 10000 Step2 : 250:5000 we use a circular mask with diameter
A general approach – 1/9
LINC/NIRVANA will require a broad spectrum of image restoration methods designed for different purposes:
- quick-look
methods, to be routinely used for a first inspection of the observed astronomical object; possibly very efficient even if not very accurate; -
ad hoc
methods for the restoration of objecs with specific features (high dynamic range, low photon counting, edges,..) It should be interesting to have a general approach for the minimization of a broad class of functionals with the constraints of nonnegativity and flux conservation (L 1 -norm). We investigate the extension of an approach proposed by Lanteri et al. (2001-2002).
A general approach – 2/9
Minimization of functionals of the following type: J(f,g) = J 0 (Af;g) + m J R (f) , where the first term is a discrepancy between detected and computed data and the second term a regularization functional, with the additional constraints: f(n) 0 , n f(n) c .
The constant c can be chosen in such a way that the flux of the restored object is consistent with the fluxes of the detected images (constant backgrounds) : c 1 p j p
1 n g j (n) b j .
A general approach – 3/9
The approach relies on the Split Gradient Method (Lanteri et al. 2002), i. e. on the following decompositions of the gradients:
f
f
f J(f;g)
U(f;g)
V(f;g) J
0
(Af;g)
U
0
(f;g)
V o (f;g) J R (f)
U R (f)
V R (f)
where U, V are nonnegative arrays (V positive). It is obvious that these decompositions always exist even if they are not unique. For practical applications their explicit expressions as functions of f and g = {g 1 , …., g p } are needed.
The approach is basically a descent method.
A general approach – 4/9
The basic algorithm is as follows: initialize with
f (
0
)
0 , satisfying the flux constraint , given
f f ( k
1
) ( k)
, compute renormaliz
f
:
( k)
a
f
e
f U
f ( k) ;g V
ω f
V ( k) ;g
ω f ( k) ;g
ω ( k
1
)
in such a way that it satisfies ; the flux constraint .
a and w are relaxation parameters; a is the step in the descent direction: it can be chosen in order to assure nonnegativity and convergence.
A general approach – 5/9
In the particular case of step 1, the algorithm takes a very simple form:
f ( k
1
)
f ( k)
U V
f f ( k) ;g ( k) ;g
ω
; convergence is not guaranteed (even if it has been proved in particular cases) but nonnegativity is automatic !
By taking into account the structure of the discrepancy functional in the multi-images case, one gets:
ω f ( k
1
)
f ( k)
j p
1
U
0
,j
f j p
1
V
0
,j
f ( k) ;g j
( k) ;g j
μ U R μ V R
.
A general approach – 6/9
The dependence of the algorithm on the multi-images suggests the following general OS-version: initialize given with
f ( k)
, set
f j (
0
)
k
0 , satisfying 1 mod the flux constraint
p
and compute : ,
f ( k
1
)
renormaliz
f ( k)
U
0
,j V
0
,j
e
f
f f ( k) ;g j ( k) ;g j
μ p μ p U V R R
ω
;
( k
1
)
in such a way that it satisfies the flux constraint .
In such a way each iteration has the same computational cost of a single-image iteration.
A general approach – 7/9
Two examples: - Poisson noise:
U
0
,j
j
A T j A j f g j
b j , V
0
,j
j
1 In the case m =0 one obtains the OS-EM algorithm.
- Gauss noise:
U
0
,j
j
A j T g j , V
0
,j
j
A j T A j f
b j
In the case m =0 one obtains an OS version of ISRA (Iterative Space Reconstruction Algorithm) with the addition of background contribution.
A general approach – 8/9
Possible regularizations:
J R
1 2
||f
f (R) ||
2
, U R J R
2
f||
2
, U R
f (R) , V R
D
D T
f, V R
f ;
I
D T D
f ;
D is a matrix with nonnegative entries. Example: the discrete Laplacian D = 4(I-D)
J R
n f
ln
f f (R)
, U R
ln
f c , V R
1 ln
f (R) c
A general approach – 9/9
As an example we obtain an acceleration of the OSEM (Ordered Subsets Expectation Maximization) method (Hudson and Larkin, 1994) derived from Tanaka, 1982: initialize given with f (0) 0, satisfying the flux constraint
f (k)
, set
j
k
1 mod
p
and compute : ,
f (k
1
)
f (k) A j T A j f g j (k)
b j ω
; renormaliz e
f (k
1
)
in such a way that it satisfies the flux constraint .
Computational cost per iteration reduced by a factor 4/(3p+1); with w =2, reduction of the number of iterations by a factor 2.
Objects with high dynamic range – 1/5
A model of young binary star (simulation based on observations of the GG Tau system)
Objects with high dynamic range – 2/5
: One possible approach can be provided by an ‘‘adaptive’’ penalization of the I-divergence (Geman and Mac Clure, 1987) J R ( f ) 2 2
n
2 f 2 (n) f 2 (n) .
Penalization is obtained by means of a semi-convex functional (the Hessian is bounded from below); the result is a ‘‘Tikhonov regularization’’ in regions where f(n) is small with respect to and ‘‘no regularization’’ in regions where f(n) is large with respect to ; in these regions the EM algorithm is dominating. Uniqueness of the minimum is not guaranteed. The previous iterative algorithm becomes very simple.
Objects with high dynamic range – 3/5
In the case of single image the algorithm is as follows initialize given with
f (
0
)
0 , satisfying the flux constraint ;
f ( k)
set
j
k
1 mod
p
and compute :
c ( k)
[(f ( k) )
2
f ( k)
η
2
]
2
f ( k
1
)
1
f ( k)
m
p c ( k)
renormaliz
A j T A j f g ( k)
b
; e
f ( k
1
)
in such a way that it satisfies the flux constraint .
Objects with high dynamic range – 4/5
Iterates of the OSEM method (linear scale) Without Regularization K=100 K finale With Regularization
Without Reg.
Objects with high dynamic range – 5/5
Ideal PSFs SR=70% SR=26% With Reg.
Concluding remarks
AIRY , version 3.0: http://dirac.disi.unige.it
A sample of open problems: - Specific methods for the reconstruction of faint objects; in these cases the photon noise may be comparable to the (Gaussian) read-out noise so that both Poisson and Gauss statistics must be taken into account.
- Development of wavelet based methods taking into account the different angular scales of astronomical objects.
- Development of powerful acceleration methods, for instance along the lines of the Biggs-Andrews approach (MATLAB).
- Domain decomposition methods for large scale images or space-variant PSFs.
References
LBT site: http://lbtwww.arcetri.astro.it
[1] M. Bertero and P. Boccacci, 1998, Introduction to Inverse Problems in Imaging, IOP, Bristol [2] M. Bertero, and P. Boccacci, 2000, Application of the OS EM method to the restoration of LBT images, Astron. Astrophys. Suppl. Ser.,
144
, 181-186 [3] M. Bertero, and P. Boccacci, 2000, Image restoration methods for the Large Binocular Telescope, Astron. Astrophys. Suppl. Ser.,
147
, 323-332 [4] S. Correia, M. Carbillet, P. Boccacci, M. Bertero, and L. Fini, 2002, Restoration of interferometric images –I. The software package AIRY, Astron. Astrophys.,
387
, 733-743 [5] M. Carbillet, S. Correia, P. Boccacci, and M. Bertero, 2002, Restoration of interferometric images –II. The case-study of the Large Binocular Telescope, Astron. Astrophys.,
387
, 744-757
[6] D. L. Snyder, A. M. Hammoud, and R. L. White, 1993, Image recovery from data acquired with a charged-coupled -device camera, 1993, J. Opt. Soc. Am., A-10, 1014-1023 [7] W. H. Richardson, 1972, Bayesian-based iterative method of image restoration, J. Opt. Soc. Am., 62, 55-59 [8] L. Lucy, 1974, An iterative technique for the rectification of observed distribution, Astron. J., 79, 745-754 [9] L. A. Shepp and Y. Vardi, 1982, Maximum likelihood reconstruction for emission tomography, IEEE Trans. Med. Imaging, 1, 113-122 [10] A. R. De Pierro, 1987, On the convergence of the iterative image space reconstruction algorithm for volume ECT, IEEE Trans. Med. Imaging, 6, 124-125 [11] H. Lanteri, M. Roche, and C. Aime, 2002, Penalized maximum likelihood image restoration with positivity constraints: multiplicative algorithms, Inverse Problems, 18, 1397-1419