Imaging and Deconvolution David J' Wilner HarvardSmithsonian CfA - PowerPoint PPT Presentation

1 / 57
About This Presentation
Title:

Imaging and Deconvolution David J' Wilner HarvardSmithsonian CfA

Description:

Imaging and Deconvolution David J' Wilner HarvardSmithsonian CfA – PowerPoint PPT presentation

Number of Views:58
Avg rating:3.0/5.0
Slides: 58
Provided by: mmschool20
Category:

less

Transcript and Presenter's Notes

Title: Imaging and Deconvolution David J' Wilner HarvardSmithsonian CfA


1
Imaging and DeconvolutionDavid J.
Wilner(Harvard-Smithsonian CfA)
  • HIA (Sub)Millimeter Observing Techniques Summer
    School
  • August 14-17, 2006

2
References
  • Thompson, A.R., Moran, J.M. Swensen, G.W. 1986,
    Interferometry and Synthesis in Radio Astronomy
    (New York Wiley) also 2nd edition
  • NRAO Summer School proceedings
  • http//www.aoc.nrao.edu/events/synthesis/
  • Perley, R.A., Schwab, F.R. Bridle, A.H., eds.
    1989, ASP Conf. Series 6, Synthesis Imaging in
    Radio Astronomy (San Francisco ASP)
  • Chapter 6 Imaging (Sramek Schwab), Chapter 8
    Deconvolution (Cornwell)
  • Lecture by T. Cornwell 2002, Imaging and
    Deconvolution
  • Lecture by S. Bhatnagar 2004, Imaging and
    Deconvolution
  • IRAM Summer School proceedings
  • http//www.iram.fr/IRAMFR/IS/archive.html
  • Guilloteau, S., ed. 2000, IRAM Millimeter
    Interferometry Summer School
  • Chapter 13 Imaging Principles, Chapter 16
    Imaging in Practice (Guilloteau)
  • Lecture by J. Pety 2004 Imaging, Deconvolution
    Image Analysis
  • Talks by many practitioners (Chandler, Wright,
    Welch, Downes, )

3
Visibility and Sky Brightness
  • from the van Citttert-Zernike theorem (TMS
    Appendix 3.1)
  • for small fields of view

    the complex visibility,V(u,v)

    is the 2D Fourier transform of

    the brightness on the sky,T(x,y)
  • u,v (wavelengths) are spatial frequencies in
    E-W
    and N-W directions, i.e. the baseline lengths
  • x,y (rad) are angles in tangent plane relative to
    a
    reference position in the E-W and N-S directions

y
x
T(x,y)
4
The Fourier Transform
  • Fourier theory states that any signal (here
    images) can be expressed as
    a sum of sinusoids
  • (x,y) plane and (u,v) plane are conjugate
  • T(x,y)
    V(u,v) FTT(x,y)
  • in this example a single Fourier component
    encodes all
  • the spatial frequency period of the wave
  • the magnitude contrast
  • the phase (not shown) shift of wave with
    respect to origin
  • Fourier Transform image contains all information
    of original image

5
The Fourier Domain
  • acquire comfort with the Fourier domain
  • In older texts, functions and their Fourier
    transforms occupy upper and lower domains, as if
    functions circulated at ground level and their
    transforms in the underworld (Bracewell 1965)
  • a few properties of the Fourier transform
  • scaling
  • shifting
  • convolution/multiplication
  • sampling theorem

6
Some 2D Fourier Transform Pairs
  • T(x,y) AmpV(u,v)

Gaussian Gaussian ? Function
Constant narrow features transform to wide
features (and vice-versa)
7
More 2D Fourier Transform Pairs
  • T(x,y) AmpV(u,v)

Disk Bessel sharp edges
result in many high spatial frequencies Ell.
Gaussian Ell. Gaussian orientations are
orthogonal in the (x,y) and (u,v) planes
8
2D Fourier Transform Pairs
  • T(x,y) AmpV(u,v)

structure on many scales
  • T(x,y) is real, but V(u,v) is complex (in
    general)
  • Real and Imaginary
  • Amplitude and Phase
  • Amplitude tells how much of a certain frequency
    component, Phase tells where
  • V(-u,-v) V(u,v) where is complex
    conjugation (Hermitian)
  • V(u0,v0) ? integral of T(x,y)dxdy total flux

9
Visibility and Sky Brightness
10
Visibility and Sky Brightness
11
Aperture Synthesis
  • sample V(u,v) at enough points to synthesis the
    equivalent large aperture of size (umax,vmax)
  • 1 pair of telescopes ? 1 (u,v) sample at a time
  • N telescopes ? number of samples N(N-1)/2
  • reconfigure physical layout of N telescopes for
    more
  • fill in (u,v) plane by making use of Earth
    rotation (Sir Martin Ryle,
    1974 Nobel Prize in Physics)

3 configurations of 8 SMA antennas at 345 GHz
12
Aperture Synthesis Telescopes
CARMA (OVROBIMA)
(e)VLA
13
Imaging (u,v) plane Sampling
  • in aperture synthesis, V(u,v) samples are limited
    by number of telescopes, and Earth-sky geometry
  • high spatial frequencies
  • maximum angular resolution
  • low spatial frequencies
  • extended structures invisible
  • irregular within high/low limits
  • sampling theorem violated
  • information lost

14
Formal Description
  • sample Fourier domain at discrete points
  • the inverse Fourier transform is
  • the convolution theorem tells us
  • where (the
    point spread function)
  • Fourier transform of sampled visibilities yields
    the true sky brightness convolved with the point
    spread function
  • (the dirty image is the true image convolved
    with the dirty beam)

15
Dirty Beam and Dirty Image
b(x,y) (dirty beam)
B(u,v)
TD(x,y) (dirty image)
T(x,y)
16
How to analyze interferometer data?
  • uv plane analysis
  • best for simple sources, e.g. point sources,
    disks
  • image plane analysis
  • Fourier transform V(u,v) samples to image plane,
    get TD(x,y)
  • but difficult to do science on dirty image
  • deconvolve b(x,y) from TD(x,y) to determine
    (model of) T(x,y)

visibilities ? dirty image ? sky
brightness
Fourier transform deconvolve
17
Details of the Dirty Image
  • Fourier Transform
  • Fast Fourier Transform (FFT) much faster than
    simple Fourier summation, O(NlogN) for 2N x 2N
    image
  • FFT requires data on regularly spaced grid
  • aperture synthesis observations not on a regular
    grid
  • Gridding is used to resample V(u,v) for FFT
  • customary to use a convolution technique
  • visibilities are noisy samples of a smooth
    function
  • nearby visibilities not independent
  • use special (Spheroidal) functions with nice
    properties
  • fall off quickly in (u,v) plane (not too much
    smoothing)
  • fall off quickly in image plane (avoid aliasing)

18
Primary Beam
T(x,y)
  • A telescope does not have uniform response across
    the entire sky
  • main lobe approximately Gaussian, fwhm 1.2?/D
    (where D is ant diameter) primary beam
  • limited field of view
  • sidelobes, error beam (usually not important)
  • primary beam response modifies sky brightness
    T(x,y) ? A(x,y)T(x,y)
  • correct with division by A(x,y) in image plane

A(x,y)
SMA 345 GHz
ALMA 690 GHz
T(x,y) large A(x,y)
small A(x,y)
19
Pixel Size and Image Size
  • pixel size
  • should satisfy sampling theorem for the longest
    baselines, ?x lt 1/2 umax , ?y lt 1/2 vmax
  • in practice, 3 to 5 pixels across the main lobe
    of the dirty beam (to aid deconvolution)
  • image size
  • natural resolution in (u,v) plane samples
    FTA(x,y), implies image size 2x primary beam
  • if there are bright sources in the sidelobes of
    A(x,y), then they will be aliased into the image
    (need to make a larger image)

20
Dirty Beam Shape and N Antennas
  • 2 Antennas

21
Dirty Beam Shape and N Antennas
  • 3 Antennas

22
Dirty Beam Shape and N Antennas
  • 4 Antennas

23
Dirty Beam Shape and N Antennas
  • 5 Antennas

24
Dirty Beam Shape and N Antennas
  • 6 Antennas

25
Dirty Beam Shape and N Antennas
  • 7 Antennas

26
Dirty Beam Shape and N Antennas
  • 8 Antennas

27
Dirty Beam Shape and Super Synthesis
  • 8 Antennas x 2 samples

28
Dirty Beam Shape and Super Synthesis
  • 8 Antennas x 6 samples

29
Dirty Beam Shape and Super Synthesis
  • 8 Antennas x 30 samples

30
Dirty Beam Shape and Super Synthesis
  • 8 Antennas x 107 samples

31
Dirty Beam Shape and Weighting
  • introduce weighting function W(u,v)
  • W modifies sidelobes of dirty beam
  • W is also gridded for FFT
  • Natural weighting
  • W(u,v) 1/?2(u,v) at points with data and
    zero elsewhere where ?2(u,v) is the noise
    variance of the (u,v) sample
  • maximizes point source sensitivity
    (lowest rms in image)
  • gives more weight to shorter baselines (larger
    spatial scales), degrades resolution

32
Dirty Beam Shape and Weighting
  • Uniform weighting
  • W(u,v) is inversely proportional to local density
    of (u,v) points, so sum of weights in a (u,v)
    cell is a constant (or zero)
  • fills (u,v) plane more uniformly, so
    (outer) sidelobes are lower
  • gives more weight to long baselines and therefore
    higher angular resolution
  • degrades point source sensitivity
    (higher rms in image)
  • can be trouble with sparse sampling (cells
    with few data points have same weight as cells
    with many data points)

33
Dirty Beam Shape and Weighting
  • Robust (Briggs) weighting
  • variant of uniform that avoids giving too much
    weight to cell with low natural weight
  • implementations differ, e.g. SN is natural weight
    of a cell, St is a threshold
  • large threshold? natural weighting
  • small threshold ? uniform weighting
  • parameter allows continuous variation between
    optimal angular resolution and optimal point
    source sensitivity

34
Dirty Beam Shape and Weighting
  • Tapering
  • apodize the (u,v) sampling by a Gaussian
  • t tapering parameter (in k? arcsec)
  • like smoothing in the image plane (convolution by
    a Gaussian)
  • gives more weight to shorter baselines, degrades
    angular resolution
  • degrades point source sensitivity but can improve
    sensitivity to extended structure
  • could use an elliptical Gaussian
  • limits to usefulness

35
Weighting and Tapering Noise
Natural 1.7x1.4 ?1.0
Robust 1.0x0.8 ?1.28
Uniform 0.9x0.7 ?1.58
Robust Taper 1.7x1.4 ?1.30
36
Weighting and Tapering Summary
  • imaging parameters provide a lot of freedom

37
Deconvolution
  • difficult to do science on dirty image
  • deconvolve b(x,y) from TD(x,y) to recover T(x,y)
  • information is missing, so be careful!

dirty image
CLEAN image
38
Deconvolution Philosophy
  • to keep you awake at night
  • ? an infinite number of T(x,y) compatible with
    sampled V(u,v), i.e. invisible distributions
    R(x,y) where b(x,y) ? R(x,y) 0
  • no data beyond umax,vmax ? unresolved structure
  • no data within umin,vmin ? limit on largest size
    scale
  • holes between umin,vmin and umax,vmax ? sidelobes
  • noise ? undetected/corrupted structure in T(x,y)
  • no unique prescription for extracting optimum
    estimate of true sky brightness from visibility
    data
  • deconvolution
  • uses non-linear techniques effectively
    interpolate/extrapolate samples of V(u,v) into
    unsampled regions of the (u,v) plane
  • aims to find a sensible model of T(x,y)
    compatible with data
  • requires a priori assumptions about T(x,y)

39
Deconvolution Algorithms
  • most common algorithms in radio astronomy
  • CLEAN (Högbom 1974)
  • a priori assumption T(x,y) is a collection of
    point sources
  • variants for computational efficiency, extended
    structure
  • Maximum Entropy (Gull and Skilling 1983)
  • a priori assumption T(x,y) is smooth and
    positive
  • vast literature about the deep meaning of entropy
    (Bayesian)
  • hybrid approaches of these can be effective
  • deconvolution requires knowledge of beam shape
    and image noise properties (usually OK for
    aperture synthesis)
  • atmospheric seeing can modify effective beam
    shape
  • deconvolution process can modify image noise
    properties

40
Basic CLEAN Algorithm
  • Initialize
  • a residual map to the dirty map
  • a Clean component list to empty
  • Identify strongest feature in residual map as a
    point source
  • Add a fraction g (the loop gain) of this point
    source to the clean component list
  • Subtract the fraction g times b(x,y) from
    residual map
  • If stopping criteria not reached, goto step 2 (an
    interation)
  • Convolve Clean component (cc) list by an estimate
    of the main lobe of the dirty beam (the Clean
    beam) and add residual map to make the final
    restored image

b(x,y)
TD(x,y)
41
Basic CLEAN Algorithm (cont)
  • stopping criteria
  • residual map max lt multiple of rms (when noise
    limited)
  • residual map max lt fraction of dirty map max
    (dynamic range limited)
  • max number of clean components reached
  • loop gain good results for g 0.1 to 0.3
  • easy to include a priori information about where
    to search for clean components (clean boxes)
  • very useful but potentially dangerous!
  • Schwarz (1978) CLEAN is equivalent to a least
    squares fit of sinusoids (in the absense of noise)

42
CLEAN
CLEAN model
TD(x,y)
restored image
residual map
43
CLEAN with Box
CLEAN model
TD(x,y)
restored image
residual map
44
CLEAN with Poor Choice of Box
CLEAN model
TD(x,y)
restored image
residual map
45
CLEAN Variants
  • Clark CLEAN
  • aims at faster speed for large images
  • Högbom-like minor cycle w/ truncated dirty
    beam, subset of largest residuals
  • in major cycle, ccs are FFTd and subtracted
    from the FFT of the residual image from the
    previous major cycle
  • Cotton-Schwab CLEAN (MX)
  • in major cycle, ccs are FFTd and subtracted
    from ungridded visibilities
  • more accurate but slower (gridding steps
    repeated)
  • Steer, Dewdny, Ito (SDI) CLEAN
  • aims to supress CLEAN stripes in smooth,
    extended emission
  • in minor cycles, any point in the residual map
    greater than a fraction (lt1) of the maximum is
    taken as a cc
  • Multi-Resolution CLEAN
  • aims to account for coupling between pixels by
    extended structure
  • independently CLEAN a smooth map and a difference
    map, fewer ccs

46
Restored Images
  • CLEAN beam size
  • natural choice is to fit the central peak of the
    dirty beam with elliptical Gaussian
  • unit of deconvolved map is Jy per CLEAN beam area
    ( intensity, can convert to
    brightness temperature)
  • minimize unit problems when adding dirty map
    residuals
  • modest super resolution often OK, but be careful
  • restored image does not fit the visibility data
  • photometry should be done with caution
  • CLEAN does not conserve flux (extrapolates)
  • extended structure missed, attenuated, distorted
  • phase errors (e.g. seeing) can spread signal
    around

47
Noise in Images
  • point source sensitivity straightforward
  • telescope area, bandwidth, integration time,
    weighting
  • in image, modify noise by primary beam response
  • extended source sensitivity problematic
  • not quite right to divide noise by ?n beams
    covered by source smoothing tapering,
    omitting data ? lower limit
  • always missing flux at some spatial scale
  • be careful with low signal-to-noise images
  • if position known, 3? OK for point source
    detection
  • if position unknown, then 5? required (flux
    biased by 1?)
  • if lt 6?, cannot measure the source size (require
    3? difference between long and short
    baselines)
  • spectral lines may have unknown position,
    velocity, width

48
Maximum Entropy Algorithm
  • Maximize a measure of smoothness (the entropy)
  • subject to the constraints
  • M is the default image
  • fast (NlogN) non-linear optimization solver due
    to Cornwell and Evans (1983)
  • optional convolve with Gaussian beam and add
    residual map to make map

b(x,y)
TD(x,y)
49
Maximum Entropy Algorithm (cont)
  • easy to include a priori information with default
    image
  • flat default best only if nothing known (or
    nothing observed!)
  • straightforward to generalize ?2 to combine
    different observations/telescopes and obtain
    optimal image
  • many measures of entropy available
  • replace log with cosh ? emptiness (does not
    enforce positivity)
  • less robust and harder to drive than CLEAN
  • works well on smooth, extended emission
  • trouble with point source sidelobes
  • no noise estimate possible from image

50
Maximum Entropy
MAXEN model
TD(x,y)
restored image
residual map
51
Example Dust around Vega
  • tune resolution and sensitivity to suit science
  • Wilner et al. 2002, ApJ, 569, L115

Vega 850 ?m Holland et al. 1998
fwhm
2.8 x 2.1 arcsec stellar photosphere
5.3x4.6 arcsec and dust blobs
SCUBA 14 arcsec
52
Missing Low Spatial Frequencies (I)
  • Large Single Telescope
  • make an image by scanning across the sky
  • all Fourier components from 0 to D sampled, where
    D is the telescope diameter (weighting depends on
    illumination)
  • Fourier transform single dish map T(x,y) ?
    A(x,y), then divide by a(x,y) FTA(x,y), to
    estimate V(u,v)
  • choose D large enough to overlap interferometer
    samples of V(u,v) and avoid using data where
    a(x,y) becomes small

density of uv points
(u,v)
53
Missing Low Spatial Frequencies (II)
  • separate array of smaller telescopes
  • use smaller telescopes observe short baselines
    not accessible to larger telescopes
  • shortest baselines from larger telescopes total
    power maps

ALMA with ACA 64 x 12 m 12 to 14000 m 12
x 7 m fills 7 to 12 4 x 12 m fills 0
to 7
54
Missing Low Spatial Frequencies (III)
  • mosaic with a homogeneous array
  • recover a range of spatial frequencies around the
    nominal baseline b using knowledge of A(x,y)
    (Ekers and Rots 1979) (and get shortest baselines
    from total power maps)
  • V(u,v) is linear combination of baselines from
    b-D to bD
  • depends on pointing direction (xo,yo) as well as
    (u,v)
  • Fourier transform with respect to pointing
    direction (xo,yo)

(u,v)
55
Self Calibration
  • a priori calibration not perfect
  • interpolated from different time, different sky
    direction from source
  • basic idea of self calibration
  • correct for antenna-based errors together with
    imaging
  • works because
  • at each time, measure N complex and N(N-1)/2
    visibilities
  • source structure represented by small number of
    parameters
  • highly overconstrained problem if N large and
    source simple
  • in practice, an iterative, non-linear relaxation
    process
  • assume initial model ? solve for time dependent
    gains ? form new sky model from corrected data
    using e.g. CLEAN ? solve for new gains
  • requires sufficient signal-to-noise ratio for
    each solution interval
  • loses absolute position information
  • dangerous with small N, complex source, low
    signal-to-noise

56
Summary
  • Calibrated Visibilities
  • ? Fourier Transform
  • Dirty Beam and Dirty Image
  • ? Deconvolution
  • Clean Beam, Sky Model, restored Image
  • ? Analysis
  • Physical Information on Source

57
Concluding Remarks
  • interferometry samples visibilities that are
    related to a sky brightness image by the Fourier
    transform
  • deconvolution corrects for incomplete sampling
  • remember there are usually an infinite number of
    images compatible with the sampled visibilities
  • astronomer must use judgement in imaging process
  • imaging is fun (compared to calibration)
  • many issues not covered today (see References)
Write a Comment
User Comments (0)
About PowerShow.com